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Summary 


By means of the depth-indicator established on a basis of pressure- 
temperature equilibrium relationship of ferromagnetic minerals, the 
depths at which plutonic rocks have been formed are estimated in the 
Japanese Islands and the surrounding area. 

Upheavals of landmasses derived therefrom show that the amount 
of the upheaval reaches more than 20km in some places in Japan and 
is very small in the Asian Continent. 

Comparison of the data of gravity anomalies, geological structure 


and upheavals makes it possible to discuss the mountain-building 
movement and its mechanism in these regions. 


1. Introduction 


There are two streams in the recent progress on rock-magnetism in Europe 
and in Australia: one is the evaluation of continental mass movement invoked by 
Clegg & others (1956, 1957) and the other, the consolidation of the polar wander- 
ing hypothesis suggested by Creer (1957) and Runcorn & others (1956). As 
suggested by these authors, the comparison of the remanent magnetization of 
rocks from North America and Europe makes it possible to connect these two 
subjects with each other. According to their interpretation, a westward drift of 
the American continent relative to Europe commenced in Precambrian times and 
the pole position of the Earth was in the Pacific Ocean in Palaeozoic times (Creer 
& others 1958; Irving 1958a). 

During the time in which these advances have been made in Britain and 
Australia, our group has made another kind of study on rock-magnetism from 
which the mountain-building mechanism can be clarified. It was found that our 
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studies bear such a close relationship to the above-mentioned results that the 


writer would like to offer several suggestions which seem to relate these 
ideas. 


2. Estimation of depth at which plutonic rocks were formed 


Recently the author (1957) studied the effect of pressure and temperature upon 
the formation of ferromagnetic mineral solid-solutions contained in rocks (mag- 
netite-ulvéspinel and hematite-ilmenite solid-solutions) and proposed the P— T 
equilibrium relationship of these minerals. His results are shown in Figure 1 (a) 
and (b). As illustrated in the diagrams, the 'Ti-concentration (or Curie point*) 
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Fic. 1. 
(a) P—T equilibrium diagram of (b) P—T equilibrium diagram of 
ilmenite-hematite solid-solution ulvéspinel-magnetite solid-solu- 
series. % shows the molecular tion series. % shows the mole- 
percentage of ilmenite in the cular percentage of ulvéspinel in 
solid-solution. the solid-solution. 


in these solid-solutions increases as the pressure and temperature at which the 
minerals have been formed are increased. This phenomenon was common in 
normal plutonic rocks which were classified gabbro, diorite and grano-diorite. A 
line OB in the diagram is drawn in such a way that all points on it represent the 
P—T conditions inside the lithosphere estimated by assuming the average geo- 
thermal gradient and the gravitational load of the rock column in the upper part 
of the Earth’s crust. When we draw the section of these diagrams along the line 
OB, the relation of the Ti-contents (or Curie point) vs depth from the Earth’s 
surface (or equivalent P—T condition) is to be that shown in Figure 2. There- 
fore, if enough care is paid to choose suitable rocks in order to get an accurate 
result of determining depth and when either Curie point or chemical composition 
of the above-said minerals in a plutonic rock is determined, the curve in Figure 2 
will indicate the depth at which rock was initially formed. 

We have sampled a number of plutonic rocks from various localities in Japan. 
By using the above-mentioned Curie point method, depth-estimation has been 
carried out at various localities in Japan and the surrounding area. The results 
so far obtained are shown in Table 1 and Figure 3. 


* Ti-contents in these minerals can be decided by measuring their Curie point. 





Mountain-building movement in Japan and its vicinity 
Table 1 


Estimated Estimated 
Locality Depth Locality Depth 


km km 
Pokchinsan (Korea) 1°2 Tongre (Korea) 3 
Ogi (Saga) 23 Yamaga (Kumamoto) 21 
Osada (Shimane) 17 Nagahama (Ehime) 23 
Unomachi (Ehime) 23 Besshi (Ehime) 23 
Shodoshima (Okayama) 22 Awajishima (Hyogo) 21 
Kibi (Okayama) 24 Yakuno (Kyoto) 24 
Maizuru (Kyoto) 23 Sakurai (Nara) 23 
Kasagi (Nara) 23 Okitsu (Mie) 24 
Toba (Mie) 22 'Tomioka (Shizuoka) 25 
Kofu (Yamanashi) 25 Funabiki (Fukushima) 
Arizitsu (Fukushima) 12 Gorintoge (Iwate) 7 
Omoto (Iwate) 7 Miyamori (Iwate) 9 
Erimo (Hokkaido) 25 Karikachi (Hokkaido) 25 
Okushibetsu (Hokkaido) 23 Nemuro (Hokkaido) 9 


3- The mountain-building mechanism of the Japanese Islands 
inferred from the depth of plutonic rocks 


Assuming that the depth of formation of plutonic rocks distributed at various 
localities be determined, we can reasonably infer how and when the Earth’s 
crust in that area was deformed. Since we sampled these plutonic rocks from 
intensely eroded and therefore from nearly horizontal land surface, the figures 
determined for the depth indicate the amount of vertical upheaval, i.e. the extent 
of erosion as it differs from place to place in this area. 
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Fic. 2.—Variation of Curie point of ferromagnetic minerals with depth. 
(a) for the series of ilmenite-hematite solid-solution. 
(b) for the series of ulvéspinel-magnetite solid-solution. 
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We can see in the figure, that the upheavals are very small in Korea but signifi- 
cantly large along the Median-line in Japan. The values become greater approach- 
ing the Japanese Islands from the west and then smaller towards the Pacific 
Ocean, resulting in maxima upheavals here and there in Japan. There maxima 
can be seen at the Kyushu Island, the Setouchi Inland-Sea area around Hiroshima, 
the area around Lake Biwa and the central mountain zone of the Japanese Alps. 


The upheavals are also great in the area of the Hidaka mountain zone of Hokkaido 
Island. 
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Since geological evidence suggests that these rocks are late Palaeozoic in age, 
the contour lines showing equal upheaval represent the topography of the crustal 
deformation which has occurred since this time. 

The writer used plutonic rocks with ferrite whose TiO2-concentration is in 
the range 5-50 per cent as test samples. When the Ti-concentration is out of 
this range the writer’s P—T diagram of titaniferous ferrites is not strictly valid, 
and his depth scale is unsuitable. 

The extent of the erosion subsequently occurring at the above-mentioned 
maxima is very great and, when averaged out, amounts to 24km. The total mass 
of land eroded and probably transported into the Pacific Ocean or the Sea of 
Japan works out to be more than 1-1 x 1016 tons. The contour lines in the figure 
show the topography of the mountains if the amount of erosion is very small. 
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Fic. 3.—Iso-upheaval map along the Japanese Islands and the 
surrounding area. 


During the Annual Meeting of Society of Geomagnetism and Geoelectricity 
of Japan held in 1958 in Tokyo, when the writer first presented his original map 
similar to Figure 3, Miyabe of Geographical Survey pointed out a marked 
resemblance between the configurations of the Bouguer gravity anomaly contour 
lines as presented by Tsuboi & others (1953/56) and those of the contour lines 
of the upheavals. Where the maxima of the upheavals appear, there exist fairly 
large negative gravity anomalies, as shown in Figure 4. This seems to indicate a 
large mass deficiency beneath the place where the mountain once existed. 

As to the negative gravity anomaly in the mountainous region, many authors 
have focussed their attention on the static model insisting on the necessity for a 
plane of compensation in the mantle. At the depth of this plane the excess gravity 
load due to the mountain would be compensated by the mass deficiency beneath. 
Some authors in interpreting the gravity anomaly found in Japan have adopted 
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this kind of compensation and have obtained the depth of the isostasy. Another 
author has suggested that the negative Bouguer anomaly in the mountainous 
region in the Middle Honshu Island is an isostatic anomaly. 

All these authors, however, in their calculations assumed the present height 
of the mountains as the cause of the above-mentioned excess mass which should 
be compensated. This height, from the present writer’s results, represents but a 
very small relic of eroded mountains, which once might have risen up to noble 
height as the Himalaya—Tibetan mountains do today. These assumptions should 
be true if isostatic balancing occurs immediately after erosion has taken place 
continuously through geological time to the present. Virtually no evidence 
supporting the above-mentioned assumption is known. 

It is more reasonable to assume that compensation occurs some time after the 
surface mass change caused by erosion because of the very viscous nature of rocks 
as has been suggested by Meinesz (the kinematic viscosity is of the order of 102° 
c.g.s.). For example in the Setouchi Inland-Sea area and Kyushu Island there is 
no high mountain at present, and the gravity has been measured nearly at sea 
level. Nevertheless, we can see a remarkable negative anomaly resulted from the 
overcompensation. 


4- Formation of the circum-Pacific mountain chain 


If the continental drift of America and Europe had occurred during Pre- 
cambrian—Mesozoic times, then the Pacific Ocean experienced great compression. 
Because the basement rocks under the Pacific Ocean and continental material are 
harder than the sediment deposited around the periphery of the continent, 
squeezing of the latter sediment in vertical directions might have been very 


severe. As a result, the circum-Pacific mountain chain could quite reasonably 
have been caused by such squeezing due to relative movement of the two conti- 
nents. Irving (1958b) has reported significant north-eastward drift of the Austra- 
lian continent, and this would have had great influence on the compression of the 
Pacific Ocean and, therefore, on the formation of the mountains. 

Thus, if the models for the formation of the circum-Pacific mountain zone 
are correct, the upheavals revealed by the Curie point measurements in plutonic 
rocks and the sinking of the sediments inferred from the Bouguer anomalies 
provide very good evidence. The squeezing of the sediments must have resulted 
in the compression of original basement materials into a sinusoidal structure. 
Such a tendency has been proved by the drilling survey near the Bikini Islands. 
A core taken from 1-5 km beneath the sea floor was found to contain a Cretaceous 
basalt and some coral relics. In the northern Pacific Area, a number of guyots have 
been found. This geological evidence illustrates the downwarp of the oceanic 
basement, whereas the Mariana folding and the Albatross ridge in the south-east 
Pacific shows the upwarping of the basements. 

The original shape of the upper bulge in Japan should have some relationship 
to the present depth of the Japanese trench. 

The tangent of the original slope towards the Sea of Japan and the Pacific 
Ocean might have been reduced by the transportation of the eroded mass from 
the continental land and its deposition on the original slope (see the schematic 
diagram Figure 5 illustrating the transfer of sediments). Thus the motive force 
of the Earth’s crustal movement occurring in Japan in the late Palaeozoic times was 
probably caused by the compression of the Pacific basement. 
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A 
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Fic. 5.—Transportation of land mass to the oceans by erosion, 
A: original land mass. B: present land surface. 
C: present slope. D: original slope. 


In Figure 6 (a) and (b), geological structures which have been inferred from a 
geological survey in Japan are shown. This interpretation has been made by Shiki 
(1960) and two sections are given, one of which shows the geology at the middle 
Permian and the other that at the late Permian. 
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Fic. 6(a).—Geological structure at the middle Permian mountain- 
building movement of the Japanese Islands (after Shiki). 
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Fic. 6(b).—Geological structure at the upper Permian mountain- 
building movement of the Japanese Islands (after Shiki). 
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The present author has been encouraged by these diagrams because they 
suggest downward compression of the geosynclinal sediments. 

If these diagrams are correct, then even without the data from the gravity 
anomalies and the evidence concerning continental drift, the compression of the 
Pacific Oceanic basement and the original mountain shape we observed become 
consistent. 
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Evaluation of the Dominion Observatory Bronze 
Pendulum Apparatus 


Angus C. Hamilton 
(Received 1961 March 10) 


Summary 

From an analysis of all the observations on the pendulum pier in 
Ottawa over a three year period it is shown that there is a long term 
trend in the periods of the bronze pendulums; there is a gradual 
increase of some 4 x 10~’s followed by a decrease of about the same 
amount which is believed to be due partly to a change in operating 
temperature and partly to creep in the bronze metal. From analysis of 
variance it is shown that the standard deviation for a gravity difference 
may be as low as 0°15 or as high as 0-58 mgal depending on the method 
of computation; this shows that the internal consistency is much better 
than the consistency between sets of observations. 

The temperature control system is discussed and the record of 
thermistor resistances for one pendulum pair is presented to show that 
there is negligible correlation between temperature and residual errors 
and only partial correlation with the long term drift of the pendulum 
periods. Other possible sources of error are discussed and it is con- 
cluded that shocks of undetermined origin affect the pendulum each time 
they are removed from the thermostated case. 


1. Introduction 


When modern gravity meters capable of making sensitive measurements pro- 
portional to gravity became available it was at first assumed that pendulum appara- 
tus would no longer be required. Although the number of pendulum observations 
needed has decreased, the need for pendulum apparatus has continued with even 
greater emphasis because of the need for an international network of control 
stations. Several research groups have undertaken development of pendulum 
apparatus with the objective of measuring gravity differences to an accuracy of 
o-1 mgal. While this objective has not yet been achieved, sufficient improvements 
have been made in several sets of apparatus that the error in measuring gravity 
differences has been reduced from one milligal down to about one-third of a 
milligal. Isolated examples suggesting much better accuracy have been mentioned 
but so far there are no published results which show beyond doubt that accuracy 
better than 0-3 to 0-4 mgal can be consistently attained. 

At the Dominion Observatory Dr L. G. D. Thompson designed apparatus 
with which to use bronze quarter-metre Mendenhall pendulums constructed 
about the turn of the century by the United States Coast and Geodetic Survey 
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Table 1 
Observations on the pendulum pier 


Room 8, Geophysical Laboratory, Dominion Observatory, Ottawa 


(In units of 10-7 s, with 4.994 400°0 subtracted from each observation) 


Pendulums Pendulums 


w | iciianiatet aenaieiieteeianaitianie 
4-5 4-6 5-6 4-5 4-6 5-6 

Date S Date S Date ‘ Date S Date Date S§ 

1956 1957 

Oct 3 62:9 Oct 5 56:7 Oct 23 59°9 Apr 18 66:4 Apr17 63:5 Apr 16 67°6 
59°5 55°0 60°3 65°8 65°9 67°8 
61°3 56°8 60°2 66°5 65°9 66°5 
60°6 67°5 65°1 69°1 

65°5 


————— 


Oct 28 61°5 56:0 
65°0 57°97 May 23 66-3. May 2669:0 May 27 67°8 
63°2 55°9 66°3 67-0 69°1 
64°1 59°6 66°5 67°0 70°4 


61°7 
62°7 
61°6 


May 30 69°6 
70°6 
68-0 


60°9 69°4 Jun 19 68:1 Jun 15 69:2 
60°6 67°3 68°8 71°4 
60°0 68°8 67°6 70°3 
61-4 67°3 
601 
66°5 67°3 71°8 
64°2 58-4 67°8 67°5 70°8 
63°9 62°4 65°7 68 +1 70°8 
63:1 66:2 
61°7 65°6 64°8 67°8 
64°7 69°3 
Mar 12 64°1 Mar 1365-0 Mar 14 67-2 66°5 
65°5 64°5 69°3 689 
62°9 63°6 68-2 69°9 
63°8 64°3 68°3 
71'0 
71°'0 
7°°9 


Mar 21 63:5 Mar 20 64:4 Mar 19 68°5 
64°1 65°8 69°2 
65°1 64°1 68°7 
65°4 63°9 69:2 69°0 

69°0 
70°'o 


69°2 
68-9 
70°2 
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4-5 
8 


65°7 


69°0 
67°7 
66°6 


Jun 26 66°8 
67°8 


July 2 659 
65:6 
666 


July 9 6571 
66°4 
65°4 


Sept 18 63°4 
62°1 
64°8 


Oct 27 62°4 
64°4 
61°3 


Nov 12 65°5 
68+1 
66:2 


Table 1—continued 





4-6 
Date § 


Apr 10 65:0 
65°5 
68-0 


69°2 
71°O 
72°4 


68°7 
72°5 
72°0 


70°2 
71°9 
70°'0 


68:8 
711 
71°4 


July 3 64:0 
67°5 
65°8 


July 10 65°6 
64°5 
65°6 


Sept 19 59°9 
60°4 
58°8 


Nov 7 63°5 
64°2 
65°2 


Nov 14 64°4 
65°0 
63°1 


ee 


5-6 
Date § 


Apr 11 68:8 
68°8 
69°2 


Apr 15 71°6 
71°6 
68-6 


Apr 16 708 
68°3 
70°2 


July 4 74°2 
70°7 
726 


July 11 74°4 
74°0 
72°9 


Sept 20 69°4 
70°'2 
68-6 


Oct 31 65°2 
67°9 
65°7 

Nov 17 70°1 


69°4 
70°6 


Pendulums 


Sean er oer re FRAME Teper ear 


4-5 
Date 


1959 

May 21 68-2 
68°5 
68:2 


Jun 3 66% 
66-4 
66:0 


July 7 63:0 
62°2 
62°4 

July 15 62:2 
62°4 
63°6 
63°2 
63°1 
62°3 


58'5 
58°9 
56°6 
62°7 
61°6 
62°4 


M = 646 


N = 89 


4-6 
Date 


May 26 67°1 
66°5 
66°2 


Jun 4 62°5 
63°6 
64°8 


July 8 62-0 
63°4 
63:6 


July 16 64°7 
63:0 
62°6 


Oct 30 63°5 
62°2 
62°6 


Nov 27 61°6 


60°0 
61°8 


N = 95 


5-6 
Date 


May 27 73°5 
72°4 
748 


Jun 5 724 
73°4 
72°38 


July 9 72°0 
714 
72'2 


July 17 70°5 
72°4 
69°6 


65'0 
65°6 
68°3 


Nov 30 64°6 
65°38 
66:2 
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Swick 1921). Tests were begun in 1956 and since then the apparatus has been 
tested repeatedly on the pier in Ottawa and during three seasons in the field 
(Thompson 1959; Winter & Valliant 1960; Winter, Valliant & Hamilton, in 
press 1961). From the start of tests to the end of 1959, no significant changes in 
procedure were introduced, the only change in equipment was the addition of a 
second set of pendulums, Numbers 1, 2, and 3, early in 1958. Although there is no 
significant difference in the performance of the two sets of pendulums, this analysis 
is limited to observations with pendulums 4, 5, and 6 used throughout the period. 
These observations provide sufficient data for an objective evaluation of the 
accuracy of the equipment and of the performance of the temperature control 
system. 


2. The observations 


The observations that have been made with pendulums 4, 5 and 6 on the pier 
in the Geophysical Laboratory at Ottawa are listed in Table 1. This includes all the 
completed observations; inevitably many observations are classed incomplete due 
to poor time signals, voltage variations, recorder failure and other difficulties. ‘The 
symbol, S, is the conventional term for the mean period of two pendulums swung 
as a pair. 

Usually three swings, each lasting about 50 minutes, were made with a pair 
of pendulums on one day, although in some of the early tests two, four or five 
swings were made on one day. At the end of a day the thermostated, evacuated 
case was opened, one pendulum of the pair was replaced by a new pendulum and 
the case was sealed, evacuated and allowed to reach thermal equilibrium over- 
night in preparation for observations the following day. 

The average of the mean periods observed on one day is designated P, thus 


P = (> S)in (1) 


in which n is the number of swings made in one day. The values of P correspond- 
ing to the data in Table 1 are plotted in Figure 1 along with histograms of their 
distribution. 

In appraising the results of these and other pendulums there is always doubt 
as to the inherent accuracy of the apparatus. Usually some type of analysis for 
internal consistency is carried out and then some comparisons are made with pre- 
viously measured and equally uncertain gravity differences. By this approach it 
is difficult, if not impossible, to arrive at an objective evaluation of the performance 
of any apparatus. On the other hand the large number of observations under 
similar conditions on one pier at various intervals over several years can be analysed 


independently and collectively to give a good picture of the performance of the 
bronze apparatus. 


3. Analysis of the observations 
(a) Long term trend 
From the plot of the points in Figure 1 it was immediately apparent that there 


was a long term trend in the periods. As a first approximation to this trend the 
data for each pair were fitted to a curve of the form 


P = aT?+bT+c (2) 
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in which P is the average period for the swings on one day, 

T is time in weeks, 

a, b, ¢ are constants to be found, 


The curves and equations are shown on Figure 1. 


bk 1956-" 1957 4 1958 —}959 ——_——-— Year 
o|NIp'3}FiMjAlmlz|a!als!oln|pia\rimlalmis|stalsiolwip|s|rimialmis{s{a!s!O{N|p [Month +—t 
0'8 £ 'M BR 0 So Tw Om ; f igo), the | 
A a to—t—p—+ 4 ft ff 4 tt tt 

















67+ 
al 
65 
64 
63 
62 ++ 
499,446) 
499, 4460 -%- ; tt ; t—+ t+ h number of residuals) 
10 
69 | 


w 
4 
x 
vv 
38 
c 
— 
is 
i 
~ 


SS SS Se ee eee ee $4 4 4 4 4 $4 tf) | 1 a 
; j | | | 
499.4459 | we j forty teghy tik ptf ptr party fb} ty 4 ° Re i 
0 8 it 2 2 7. SO 88, OF HR OO 8 Ae 415 Oe | OB Week Residuals x10""s 
OND J\F.MAMJ/J AS|OND J\F.MAMJ|J/A\$\0\NiD! Month . 
~ 1958 ~- 1959 


++ ++ 
0-4) 20 0 + +40 +00 


OINID JiFMAMJ JIA'S 
1» 195 bo. 1957 


+ Year 


Fic. 1.—Observed periods, Ottawa, 1956-59 for pendulums 4, 5 and 6. 


There is no evidence to show conclusively what is causing this long term drift 
in the pendulum periods. As discussed subsequently and illustrated in Figure 3, 
the readings of the thermistors that monitor the temperature in the pendulum 
case are changing gradually but whether this represents a change in temperature 
or a drift in thermistor characteristics is uncertain. Even if the temperature within 
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the case is changing as much as indicated by the lid thermistors there is still a drift 
of some 2 x 10-7 s per year to be explained. It is believed that most of the drift is 
due to creep in the bronze metal similar to changes that have been noted for 
length standards (Field 1939). Analysis of the standardization periods of pendu- 
lums 1, 2 and 3 for the years 1921 to 1927 (Miller 1929) shows that during those 
years there was a decrease in period averaging 6-41, 1-65 and 1°41 x 1077s per 
year respectively. A period change of 2 x 10~’s is equivalent to one part per million 
in the length of a quarter-metre pendulum. Several factors, such as rounding of 
the knife edges and sway of the support may contribute to an increase in effective 
length but creep in the metal itself is the most probable explanation for a decrease 
in effective length of the pendulums. 


(b) Transformation of the data 
The long term trend described above changes so slowly that it does not seriously 
affect relative gravity measurements made during an interval of a few weeks but 
for an analysis over the three-year period it must be taken into account. Rather 
than deal with residuals relative to the curve a transformation* can be made such 
that the points have similar residuals about their mean value. For average periods 
this is given by 
P’ = M+P-—(aT?+6T+<c) (3) 


and for mean periods it is given by 
S’ = M+S—(aT?2+bT+c) (4) 


in which a, b,c, Pz | S§ are as defined in (1) and (2) above and M is the mean, 
of all the periods with one pair of pendulums both before and after transformation, 


thus 
M = (SS)|N = (SSN (5) 


where N is the total number of swings with one pair of pendulums. 


(c) Variance analysis 

(i) Variance (0?) of a set of swings.—For standardization and field measurements, 
observations are usually made on three successive days with the three pairs of 
pendulums; this is sometimes called a set of swings. From Table 1 it is apparent 
that there were numerous exceptions to this procedure during the testing period 
in 1956 and 1957. As a check on the internal consistency of the observations the 
variance of each set of swings with the three pairs is computed (cf. Cook 1950; 
Garland 1953; Winter & Valliant 1960; and others) by the formula 


3 3 
[> > (Sei—Pe)?/( > m-3) (6) 
k=1 i=1 k=1 


in which S, P, and n are as defined above and k = 1, 2, 3 refers to pendulum 
pairs 4—5, 4-6 and 5-6. 

By this method only one variance is obtained for the periods of the three pen- 
dulum pairs in each set. The average for all complete sets is listed in row A, 
Table 2. It follows that some of the incomplete sets of test swings are not included 
in this average. 


* To illustrate the procedure an analytical method of transformation is described although a 
graphical method was actually used. 
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This variance is a useful criterion of the care with which the measurements 
were made and of the performance of the ancillary equipment but as will be shown 
later it is of negligible value as a measure of the precision of gravity measurements. 


(ii) Variance of observed periods.—From the data in Table 1 the variance of each 
pair of pendulums was found by the formula 


h n 
[> > (Sa-M?1/(N-1) (7) 
j-l i=l 
in which S, M, n and N have been defined above and fA is the number of days, 
hence the number of sub-groups of observations. 


Similarly the variance of P about the mean has been obtained from 
h 
[ 2 (Ps — M)?}/(A—1). (8) 


These variances are listed in Table 2, rows B and C respectively. ‘These values 
are vitiated by the long term effect and consequently are not valid criteria of the 
day-to-day performance of the apparatus. Comparison of B and C with D and H 
respectively shows that the variances are significantly smaller after the effect of the 
long term trend has been removed. 

(iii) Variance of transformed periods—When a large number of observations on 
one variate are divided into groups and a mean is obtained for each group the re- 
lation between the observations, the mean, and the group means is given (Weather- 
burn 1947) by 


h 


n; h n h 
> > ('n-MP= > & (Su —Pi P+ > n(Pj'-—MP (9) 
j=l it=1 j=l i=l j=l 


in which S’ is the transformed mean period 
Mis the grand mean of all periods for one pair 
P’ is the transformed average of the periods for one day, i.e. for one 
group 
h_ is the number of groups, i.e. the number of days 
n is the number of swings on one day. 


When each term in this equation is divided by the appropriate number of 
degrees of freedom three unbiased estimates of variance are obtained which 
satisfy the identity 


(N—1)o2 = (N—h)o? + (h—1)o?. (10) 


These variances are listed in Table 2, rows D, E and G. 

If the classification into groups, such as these above, is non-selective, the three 
estimates of variance will approach equality in the limit; for finite samples we can 
test for the existence of selective classification by comparing ratios of these vari- 
ances with ratios predicted by theoretical statistical methods. 

(iv) Variance ratio tests—If a large sample, say 1 000, is randomly drawn from a 
normally distributed population that is several times larger than the sample, one 
would expect the variance of the sample to differ very little from that of the popula- 
tion; for small samples from much smaller populations an appreciable difference 
can be expected in the variances hence the ratio of these variances will depend on 
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the size of the sample and of the population. Tables of variance ratios (F’) for nor- 
mal populations have been published by Snedecor (1934). For the number of 
degrees of freedom of the sample and of the population these tables list two values 
designated as 5 per cent and 1 per cent points of F. The 5 per cent point for 
instance is that value which, for the given degrees of freedom, will be exceeded in 
5 cases out of a hundred. These values can be used to test the hypotheses that 
samples are either from the same normal population or are from normal popula- 
tions of equal variance. A value of F less than the 5 per cent point is not con- 
sidered significant, a value greater than the 1 per cent point is highly significant 
‘and throws doubt on the truth of the hypothesis, hence these points are also 
known as significance levels. 

Using the transformed observations from which the long term trend has been 
removed we have, for each pendulum pair, the variance from a total population, D, 
(Table 2) the variance from groups, £, and the variance from group means, G. 
As shown by (g) and (10) these variances should be equal. By comparing their 
ratios with predicted values of F we can test whether or not the differences in 
variance are significant. 

In the lower section of Table 2 the variance ratios, degrees of freedom and 
corresponding significance levels of F are summarized. As all the ratios are much 
greater than the 1 per cent of F there is no doubt whatever that the differences in 
variances are highly significant. 

(v) Why do the variances differ so widely?—As far as physically possible every 
swing has been made under similar conditions; the same temperature, pressure, 
initial arc, etc., were maintained for all observations. Yet it is generally assumed— 
and this analysis confirms it quantitatively since og? < op?—that two swings on any 
one day tend to differ less than two swings on different days. The variance ratio 
D/E averages 5 whereas the probability is only o-o1 that it should exceed 1-8. As 
it should, variance E does not differ greatly from variance A; both are based on 
deviations from group, or daily, means, hence both are a measure of internal con- 
sistency. Variance D is based on deviations from the grand mean, hence it is a 
measure of external consistency. As shown by equations (9) and (10) variances 
D and E will be approximately equal if there is no selective classification operative. 
Obviously the grouping by days is highly selective; in other words there is a dis- 
crepancy between internal and external consistency. 

The suspected cause of this discrepancy is attributed to shocks received during 


the transfer of the pendulums from the sealed case to the heated carrying case and 
back between successive days. 


(d) Serial correlation tests 


From the variance analysis above it was apparent that a pendulum observation 
is not independent of its predecessor, in other words the observations are serially 
correlated. To confirm this indication serial correlation tests were made on 5S, P, 
S’, and P’; these tests are summarized in Table 3. 

For these tests it is necessary to compute a coefficient 


(11) 


in which 8? is the mean of the sum of the squares of first differences and o? is the 
variance. 
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For example 


h-1 
> (Psj11—P3)?)/(A—1) 
j=l 





(12) 


age 
[2 (Ps— MPA 1) 


As for variance analysis, tables giving the 5 per cent and 1 per cent points of K for 
different degrees of freedom are available (Ferber 1949).* In this case a value of K 


less than 1 per cent significance level indicates the existence of positive serial 
correlation. 


Table 3 
Serial correlation tests 


Pendulums Significance 
4-5 4-6 5-6 Level 
I per cent 5 per cent 


0°36 0:28 
0°53 0°46 
°o-6 06 
I°t °o'9 


In Table 3 significant positive serial correlation is indicated for both S and P. 
This is to be expected because of the long term trend. Even after the transforma- 
tion which removes the effect of the long term trend there is still significant corre- 
lation for S’ and marginal correlation for P’. The significant positive correlation 
for 5’ confirms the evidence from variance analysis that on any one day the obser- 
vations are clustered together. The marginal correlation for P’ can be interpreted 
either as an indication that a more complex curve is needed to describe the long 


term trend or that some of the clustering effect noted for S’ carries over from day 
to day. 


(e) Comments on the behaviour of the pendulums 

From a study of the information in Figure 1 and Tables 2, and 3, we can make 
some deductions concerning the behaviour of the bronze pendulums. 

(i) As long as the pendulums are left in the thermostated, evacuated case and 
only raised and lowered between swings the periods are quite consistent. This is 
confirmed by the low variances A and E. 

(ii) The effect of opening the case, removing the pendulums to the heated 
carrying case, and returning them to the sealed case introduces shocks of unknown 
origin which tend to increase the variance by a factor of five. 

(iii) These shocks are essentially random in size and direction but, as indicated 
by the serial correlation of P’, each one is influenced slightly by its predecessor. 

(v) There is no evidence to show that the pendulums will recover from these 
shocks while in the sealed case. There are insufficient data to test this statistically 
but when the pendulums have been left in the sealed case between sets of swings 


* Note.—In Table 3 the values of the significance levels for S and S’ have been obtained by 
extrapolation from the Table on p. 525 of Ferber. 
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(bracketed groups, Figure 1) regardless of the time lapse they appear to behave as 
we would expect if the swings had all been made as one large set. 

(v) There is also no evidence to show that long field trips or long periods of 
quiescence have any more effect than a simple transfer to the carrying case and 
back to the sealed case. This can be interpreted as a tribute to careful handling 
and constant heating on field trips. 


4. Standard deviation for a gravity difference 
(a) Observations at Ottawa 


For each average variance in Table 2, corresponding standard deviations have 
been tabulated and the standard deviation of P’ has been found on the assumption 
that there are three observations in each group, thus 


o(P’) = o(S’)/3?. (13) 


This is not strictly accurate but it is adequate for this type of analysis. The equiva- 
lent standard deviation of P’ in terms of milligals listed in the second last column 
of Table 2 has been computed from the approximate relation 


o(P’) (mgal) = 0-4 o(P’) (s x 1077). (14) 


In measuring a gravity difference between two stations six pendulums are 
grouped in two sets of three so that observations are made with six pairs of pen- 
dulums. If o(P’) is the standard deviation in milligals of the observations with one 
pair of pendulums at one station then by the laws of error propagation the standard 
deviation of the difference between two stations is given by 


o(Ag) = o( P’)6-42#. (15) 


Values of o(Ag) corresponding to the different variances are listed in the last 
column of Table 2. Using the same observed data, standard deviations differing 
by a factor of four have been obtained; this illustrates the importance of defining 
clearly the method by which standard deviations for gravity differences are com- 
puted. Unfortunately it is not possible to say that one of these values is correct 
and the others are wrong; each contains useful information about the apparatus 
and studying them as a group yields valuable information about the performance 
of the apparatus. 

(i) Standard deviations from rows A and E of Table 2 are equal. This is to be 
expected as they are derived from residuals of periods S, or S’ about group means 
P and P’. The low value of 0-15 mgal shows that the deviations about the group 
means are small hence the internal consistency is good. This value of standard 
deviation would apply if a gravity difference could be measured without exposing 
the pendulums to the shocks which they apparently receive each time they are re- 
moved from the evacuated thermostated case. 

(ii) Standard deviations from rows G and H are 0-58 and 0-57 mgal respectively ; 
although both are derived from the residuals of group means, P’, about the grand 
means, M, the small difference arises because the values in row G are weighted 
proportional to the number of observations in each group whereas in row H each 
group is given equal weight. This is the order of magnitude for the standard 
deviation of a gravity difference that would apply if field observations were made 
without standardization swings before and after each field trip; in other words it is 
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the value that would apply if field procedure were not designed to utilize the 
beneficial effects of serial correlation. 

(iii) The standard deviation of 0-35 mgal given in row D is based on the residuals 
of the transformed periods, S’, about the means, M, for each pair of pendulums. 
This is the most direct and hence the best method by which to evaluate the in- 
herent accuracy of the apparatus. It is fortuitous that the mean of the values from 
(i) and (ii) above, is 0-36mgal. 


(b) Internal analysis of field observations 


In the field usually only one set of observations is made at each station hence 
only the variance analysis by the formula in row A, Table 2, can be applied. This 
has been done for each season’s work and the results summarized in Table 4. It is 
apparent that the internal consistency of the results in the field is comparable with 
that on the pier in Ottawa. 


(c) Pair analysis of field observations 


It is customary to compute the gravity difference between two stations inde- 
pendently for each pair of pendulums. The standard deviation can then be com- 
puted for a pair and also for the mean gravity difference. Averages for each season 
are summarized in Table 4. 


Table 4 


Sammary of standard deviations and comparisons of gravity differences (mgal) 


Internal Pair Comparisons Ottawa 
Location Year Analysis Analysis (without regard Pier 
o Ag oAg to sign) o Ag 


Canada and 1957 o-17* : 0°33 
(Only one set 
United States of 3 
pendulums used) 


Prairie Provinces 1958 0-16 0°20 
Canada and 

Europe 1959 0°14 0°16 
Ottawa 1956-59 ols 


* Two erratic observations omitted in computing this value. 


(d) Comparisons with other values 


Each gravity difference that was measured with these pendulums has been 
measured many times with gravity meters and in some cases with other pendulum 
apparatus. From these previous measurements a value—usually the mean—has 
been adopted in order to make a comparison with these results. Individually each 
of these previous values is of dubious accuracy but on the average they are not 
likely to be seriously in error. It is interesting to note that these comparisons, 
listed in Table 4, do not differ significantly from the standard deviation derived 
from analysis of the results on the Ottawa pier. 
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(e) Weighting observations for network adjustments 

From the discussion above and from data in Tables 2 and 4, it follows that there 
may be several distinctly different values for the standard deviation of a gravity 
difference. Thus there are several approaches to the weighting of the observations. 
Weighting frequently is based on the internal consistency of each set of measure- 
ments. This appears to be scientifically rigorous but unless it can be shown that 
the internal analysis takes account of all the factors it is virtually meaningless. 
From this investigation it appears that the internal analysis takes account of all 
small errors and omits a large one. 

Differences measured with this apparatus have a standard error of the order of 
0-35 mgal and there is no way of determining from any one set of measurements 
whether the error is greater or less than this value. If the measurements have been 
made with the usual care—o(S) of the order of 1-1 x 10-7s—there is no way of 
deriving a figure for the error in the gravity difference; whether o($) is o-g or 
1-3 x 10~’s is insignificant when it is known that the whole set may be in error by 
20x 10~’s. Before the results with different apparatus can be weighted validly 
for adjustment of the first-order world network a similar objective appraisal of each 
apparatus is needed. 






































Fic. 2.—Location of thermistors in case wall relative to pendulums. 


5- Pendulum temperature tests 


(a) Temperature control and monitoring system 

Variations in temperature inside the pendulum case are limited to 0-o1 degC 
by a thermistor (Western Electric Type 10A) embedded in the wall of the case 
acting as a thermostat. Selection of the operating temperature (about 40°C) is 
arbitrary but once selected must not be allowed to change. To check that the 
temperature remains stable and uniform within the case seven thermistors (Western 
Electric Type 12A) are mounted in it; five are on the wall and two on the lid (see 
Figure 2). The resistance of these thermistors can be measured by a bridge 
circuit to give the temperature at any sensing point at any time. The coils that 
supply the heat are wound such that in any horizontal plane the temperature is 


radially symmetrical but very sensitive adjustments are necessary to eliminate 
vertical gradients. 
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(b) Temperature sensitivity of the pendulums 
A change in temperature of o-o6degC will cause a change in period of 2-5 x 
10~’s, equivalent to an apparent change in gravity of 1mgal (McDiarmid 1915). 
As temperature increases, length of the pendulums increases, periods increase and 
observed g decreases. Thermistor resistance decreases with increase in temperature ; 
for these thermistors the coefficient is very nearly —o-oosdegC/ohm. ‘Thus for a 
change in resistance of 1ohm: change in period is 


A I X —0°005 ~. 
p = ————- = -0°2x 107’s (16) 
0°024 


and change in gravity is 
IX —0°005 
— 0°06 


(c) Thermistor resistance for pendulum 5-6, Ottawa 1956-59. 

In Figure 3 for each thermistor the average resistance recorded during swings 
with pendulum pair 5-6 at Ottawa is plotted. Fluctuations in resistance recorded 
during 1956 and early 1957 reflect adjustments that were made in the heat control 
circuitry to eliminate vertical gradients. Subsequent fluctuations in the record 
that correspond to case overheat or to a breakdown of the control circuitry are 
noted on the figure. 


Ag = 0-08 mgal. (17) 


(d) Average deviations in thermistor readings 

For the period of this analysis the average resistance for each thermistor was 
found and from these averages were found the individual deviations at each point 
and the average deviation of the case temperature from normal. For this purpose 
the thermistors were grouped three ways: 


(i) All seven thermistors 1, 2, 3, 4, 5, 16, 17, 
(ii) ‘Two lid thermistors 16, 17, 
(iii) Five case thermistors, 1, 2, 3, 4, 5, 


the data from each group are plotted in Figure 3. An analysis to detect a long 
term trend was made in the same manner as in ‘““Long Term Trend” above. The 
curves and equations are shown on the Figure. It is apparent that, to a limited 
extent, the trend corresponds with the long term trend in pendulum periods. 


(e) Correlation between thermistor resistance and pendulum periods 
The temperature deviations for each of the three groupings above were tested 
for cross correlation with the deviations of pendulum periods (Figure 1). From 
the temperature sensitivity of the pendulums it might be expected that a large 
positive (negative) AP would correspond to a large negative (positive) AR. This 
was not confirmed. Cross correlation was tested by the formula 
N 1 N N 
> APnARn—— > APn > ARa 
n=1 N n=l n=1 
Stas: * a N 1 (18) 
[{ > APn?—— ( > APa)*} {> ARn?——( 5 AR,)*}} 
N n-1 n=1 N ‘n=1 


n=l 





in which AP is the deviation in pendulum period 
AR is the deviation in thermistor resistance 
N = 29. 
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For the three thermistor groupings we get 
Tr, = 0°43 (1 to 17) 


ry; = 0°24 (16 & 17) 
rij = 0°32 (1 to 5) 


Although all three coefficients indicate negligible correlation it is interesting to note 
that the value from the lid thermistors is even lower than the others. 
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Fic. 3.—Thermistor resistance readings recorded during swings of 
pendulum pair 5 and 6, 1956-59. 


(f) Assessment of the temperature control system 

The evidence from this analysis indicates that the pendulum periods are 
affected less by temperature changes than would be predicted using the coefficient 
for temperature sensitivity. The long term trend in temperature is in the proper 
direction but not large enough to account for the trend in pendulum periods. 
There is still the possibility that the thermistors—including the control thermistor 
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—are drifting in the same direction at the same rate; this seems highly improbable 
but it could effectively mask a change in operating temperature. An absolute 
calibration of the thermistors was made at the end of 1959 but due to uncertainty 
in identification it cannot be correlated with the absolute calibration made in 1956. 
From the data for all the thermistors over the entire period the standard 
deviation of thermistor resistance is 1-°80hm; combining this with the sensitivity 
factor of o-o8mgal/ohm gives a standard deviation due to temperature changes of 
o-14mgal. This figure is high since it includes the fluctuations in lid temperature as 
well as those due to adjustment in 1956-57. Using the data from the case ther- 
mistors only, the standard deviation due to temperature changes for a set of 
measurements with one pair of pendulums becomes 0-10mgal. This appears to 
be a fair assessment of the temperature control system. 


6. Possible sources of error other than temperature 


In addition to temperature changes there are numerous factors which may 


contribute to the error of an observation. Several of these are discussed and their 
effect is estimated. 


(a) Recording 


Irregularities in the rotation of the photographic drum and shrinkage or ex- 
pansion of the photographic paper during processing can contribute at the most 


only very small errors. The maximum contribution to the standard deviation is 
estimated at 0°03 mgal. 


(b) Scaling and interpolating records 


These errors are due to small variations in the size of the dots on the record and 
to the uncertainty in estimating the centre of a dot precisely. The contribution 


to the standard deviation of a set of swings from errors in recording and scaling is 
estimated to be 0-ogmgal. 


(c) Are correction 


The initial arc is maintained at 0-0075 + 0°0003 radian for each swing and a 
small correction is applied to reduce the period to that of a pendulum swinging 
with an arc of 0-0075 radian. This correction is determined by measuring photo- 
graphically the amplitude of a reflected light beam at a distance of one metre from 
the pendulum case. A pendulum arc of 0-0075 radian gives a deflection of 30-0omm 
to the light beam. Deviations from 30-00mm are corrected for by application of 
the empirical formula C = 1-2(30-00—.x) in which x is the observed deflection of 
the light beam in millimetres and C is the correction in 10-78. The term usually 
amounts to about 0-5 x 10~? and never exceeds 1-0 x 1077; the error in deriving this 
correction should not exceed 0-1 x 10~’s for one swing or 0-1/3', equivalent to 
0:06 x 1077s, for three swings. 


(d) Knife edges and flats 

When pendulum 5 is swung on knife edge No. 1, its average period is longer 
by 56x 1077s than when it is swung on knife edge No. 2. This shows that there 
is a knife edge effect due either to rigidity or rounding. However, since it has not 
changed appreciably during three years’ observations it is fair to assume that it will 
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not have any significant effect on gravity differences. This subject is discussed in 
detail by Winter & Valliant (1960). 


(e) Timing 

For time interval measurements a portable frequency standard is aligned with 
WWY several times during a swing. It is obvious that the accuracy of the portable 
standard can be no better than the value against which it is compared. Ionospheric 
effects on the sky wave (Shaull 1950) can at times reduce the accuracy of the re- 
ceived frequency to a few parts in 10’ which can cause errors as large as 0-5 mgal 
in one swing. This large error could conceivably occur on three successive swings 
but the probability of this happening is low. In fact, it is more probable that the 
error during one swing would be compensated during a successive swing. On this 
basis the standard deviation in a set of three swings is estimated to be no larger than 
o-2mgal. 


7. Summary and conclusions 


(a) The periods of the bronze pendulums are changing slowly. Part of this 
change may be due to change in temperature as a result of drift in the control 
thermistor but it is believed that at least part of the change is due to creep in the 
bronze metal. 

(b) The standard deviation derived by combining the contributions from all 
known sources of error agrees well with the standard deviation derived from 
analysis of variance for internal consistency. Using estimates derived from the 
discussion of error sources yields the following standard deviations in g for a set 
of three swings with one pair of pendulums: 


Temperature o-1romgal 
Recording 0°03 
Scaling and interpolating records 0:06 
Arc correction 0°03 
Knife edges and flats 0°05 
Timing 0°20 


from which 


n 
og = ( > a*)! = o-25mgal. (19) 
i~1 


From analysis of variance (Table 2, row A) the standard deviation for the same 
number of swings is 0-26mgal indicating that there is no serious discrepancy in 
the assessment of errors. 

(c) The standard deviation for a gravity difference measured with all the 
Dominion Observatory bronze pendulums varies from 0-15 to 0-58 mgal depending 
on the method by which it is calculated. Thus illustrates the difficulty in getting an 
objective assessment of any pendulum apparatus; an illusion that very high accuracy 
is being attained can be created if the analysis is based on internal consistency only. 
Internal analysis provides a good criterion of observational technique and of the 
performance of the ancillary equipment but it is virtually valueless as a basis for 
assessing the accuracy of results. It is suggested that a standard deviation of 
0-35 mgal be used for weighting observations made with this apparatus. 
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(d) From this analysis it appears that the pendulums are quite consistent as 
long as they are not removed from the thermostated evacuated case whereas even 
a brief transfer to the carrying case may cause a change in period large enough to 
account for the fourfold difference in standard deviation that is found. There is 
no indication of the cause of this effect; whether it is due to thermal, physical, or 
other shocks there is no way of deciding at present. When the electronic recording 
and timing system now being developed is completed it is expected that errors 
other than those due to the pendulums themselves will be negligible; it is then 


proposed to undertake a program to locate the cause of the erratic behaviour of the 
pendulums. 
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Interaction Between Tide and Surge in the Thames 


J. R. Rossiter 
(Received 1961 March 16) 


Summary 


Observations of storm surges in the Thames reveal an interaction 
effect between tide and surge which is statistically significant, and takes 
the form of an amplification of surge height (whether positive or nega- 
tive) on the rising tide. Representing the estuary by a channel of con- 
stant depth and exponentially varying breadth the distributions of 
elevation and current along the estuary have been computed for different 
combinations of input (at the seaward end) of tide and surge. The 
results agree remarkably well with the phenomenon as observed in 
nature; they provide the basis for a simple physical explanation of the 
interaction, in which it is argued that the principal mode of interaction 
is by the surge changing the phase of the tide, and vice-versa. This 
argument requires the tides in the Thames to be of progressive wave 
type; this is shown to be true, in general, for tides subject to friction 
in an exponential estuary of constant depth. 

The computations also reveal that maximum levels are propagated 
upstream so that a higher maximum is amplified less than a lower 
maximum. 

The paper confirms a number of Proudman’s earlier theoretical 
results on interaction in estuaries. 


1. Introduction 


In continuation of the research programme conducted by the Liverpool Obser- 
vatory and Tidal Institute into storm surges around the coasts of the British Isles, 
once the basic surge-forecasting formulae for North Sea stations had been provided 
to the Advisory Committee on Oceanographic and Meteorological Research of 
the Ministry of Agriculture, Fisheries and Food (Rossiter 1959), attention was 
focused upon the problem of propagation of surges up the River Thames. 

The Flood Warning Organization, responsible for the practical task of day to 
day forecasting of danger levels, had early encountered considerable difficulty in 
estimating the amplification of surges propagated from the North Sea into the 
Thames. In particular, positive surges observed near high water in the Lowestoft/ 
Harwich area and travelling southwards with the tide, rarely turned out to be as 
dangerous as might be expected; in consequence there arose the possibility of an 
excessive number of false alarms being issued. 

This problem was considered to be a particular aspect of the larger one of inter- 
action between tide and surge. The hydrodynamical equations governing the 
propagation of long waves in shallow water, containing, as they do, non-linear 
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terms, are clearly the reason for secondary oscillations being produced when tide 
(primarily a continuous oscillation of period 12h) and surge (a transient oscillation 
of approximate ‘“‘period” 24h) co-exist. Formal expressions of limited scope for 
the secondary effects in the case of continuous oscillations of one or more periods 
can be obtained (e.g. Proudman 1953, § 145), but when transients are involved the 
most satisfactory approach is by the numerical solution of type problems. 

So far as the author is aware only two names have hitherto been associated 
with this specific problem of interaction, those of Proudman (1955 a and b, 1957, 
1958) and Doodson (1956). Proudman’s papers developed formal solutions of 
great interest; most of them included bottom friction proportional to the square 
of the current; the effects of the non-linear inertia term and the variation in cross- 
sectional area of an estuary were considered. Despite mathematical complications 
which prevented a general solution being found, Proudman’s principal 1955b 
result, that for the same sequence of meteorological conditions over the sea the 
surge up an estuary is less at the time of high water than it is at the time of low 
water, is confirmed in the present investigation. Doodson’s work was primarily 
numerical in nature; he developed a method of computation and applied it to the 
case of a uniform gulf of depth 128 ft and length 105 miles. He found no appreci- 
able interaction between a semidiurnal tide and a diurnal surge, almost certainly 
because the depth chosen was too great. 

The present paper is a numerical investigation of the propagation of tide and 
surge in an estuary similar in dimensions to the River Thames, with a view to 
explaining specific phenomena observed in that river. 


2. The nature of interaction in the Thames 


From the head of the tidal Thames (Richmond Lock) to its seaward limit is 
some 68 miles; Southend is 7 miles within this limit, and tidal records for this 
station have been taken as typical of the estuary. Hourly values of surge height for 
all important surges, both positive and negative, which were experienced during 
the years 1928 to 1938 were available from a previous investigation (Rossiter 1959). 
From these were extracted the surge height at predicted high water time (H), 
predicted low water time (L), and at the two intermediate times on the rising (R) 
and falling (F) tides; care was taken to extract surge heights in multiples of 4 from 
each tidal cycle within a surge period. 

Table 1 tabulates the number of surge heights exceeding chosen values for the 
four states of tide for (a) positive surges, (b) negative surges and (c) both positive 
and negative surges. 

This table clearly reveals a tendency for surges to be greatest on the rising tide, 
irrespective of the sign of the surge; there appear to be no significant differences 
between the distributions at the three other states of tide. These results confirm 
a similar tendency noted by Doodson (1928) on the basis of many fewer data. It is 
nevertheless desirable that the uniqueness of surge distribution on the rising tide 
be proved statistically significant, and this was done by means of the ,?-test. 
Assuming that the distribution to be expected at any state of tide be represented 
by the average of the distributions for H, F, L and R, the test shows that the proba- 
bility of random errors in sampling being responsible for the observed R distri- 
bution is less than o-1 per cent for negative surges and much less than this for 
positive surges. It is therefore clear that the observed phenomenon is statistically 


significant. 
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In order to arrive at a quantitative estimate of the observed effect, values of 
surge height (residual) were listed in units of a tidal cycle and grouped according 
to the mean residual of each cycle. In Figure 1 are plotted the mean residuals for 
each of the four states of tide as ordinates and the mean residual for the complete 
tidal cycle as abscissae. Assuming linearity over the range involved, it will be seen 
that on the rising tide surges are increased by approximately 25 per cent over those 
at any other state. Consistency between positive and negative surges is main- 
tained. 





+ at low water } af rising tide 
¢ 


5 











T ¢ a 
at high wafer = af falling fide 


+ 


5 


+ 








Fic. 1.—Values of surge height at Southend (ordinates) at four states 
of tide relative to surge height averaged over the tidal cycle 
(abscissae). Units are feet. 


An example of a typical surge which exhibits the interaction effect is given in 
Figure 2. This diagram (2a) gives the Southend residuals for the surge of 2/3 
December 1929, and clearly indicates both the trough and the subsequent peak 
occurring on the rising tide. Application of a graphical smoothing process to the 
residuals produces a damped oscillation of period near to 8h (2b), the first trough 
of which (21.30h December 2) gives rise to the very sharp trough in the surge, 
and the second major crest of which (10.00h December 3) produces the sharp 
peak in the surge. This particular example is quite typical of surges at Southend 
and stations higher up the Thames, and has been chosen because it was preceded 
by quiet meteorological conditions. The synoptic situation was such that strong 
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southerly winds spread westwards into the whole of the North Sea on December 2, 
veering to southwest in the afternoon and moderating by the morning of Decem- 
ber 3. The main positive surge may be described as a “return” surge, generated by 
the return, primarily down the East Coast, of water previously expelled from the 
North Sea by the southerly winds. 

Before considering the most likely physical mechanism for the observed 
phenomenon, certain possibilities may be excluded. In 1928 Doodson stated: 
‘General reasoning indicates that the effects of an increase of depth of somewhat 
long duration over a wide area will be to decrease the range of tide; that is, the 














. 











Fic. 2. (a) Surge heights for 1929 December 2 and 3, at Southend; 
(b) damped oscillatory component of surge. H and L denote tidal high 
and low waters respectively. 


apparent meteorological effects (as shown by the residues) are a general average 
increment of level, with minima at high water and maxima at low water. These 
conditions hold when there has been a long-continued NW wind. Apart from this 
purely tidal effect, the local wind will probably have most influence at low water 
and will tend to give maximum results between then and high water. Hence there 
are two theoretical reasons for believing that there is a marked tendency for sea 
level to be most raised, by meteorological causes, shortly after low water.” 

With regard to Doodson’s first suggestion, residuals at Aberdeen, Tynemouth, 
Immingham and Lowestoft reveal no evidence of interaction to any perceptible 
degree, strongly suggesting that the Southend effect is a relatively local one not 
associated with forces acting over a wide area. Nor does the suggested mechanism 
seem quantitatively possible in this case. During the exceptional surge of 1953 the 
mean level of the North Sea as a whole was raised a maximum of approximately 
2ft (Rossiter 1954), so small a proportion of the mean depth of the sea and its 
approaches that its effect on the range of tide could hardly produce the amount of 
interaction observed at Southend. 
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The second of Doodson’s suggested mechanisms for interaction is based upon 
the well-known theory that a wind generated water gradient is inversely pro- 
portional to the depth of water. ‘Tomezac (1952a, b) has shown this effect to be of 
some importance in German North Sea estuaries. For a 25 per cent increase in 
the water gradient down the North Sea, as observed at Southend on the rising 
tide, the theory would require a 25 per cent decrease in the mean sea level due to 
the state of the tide. This effectively restricts the area in which the theory could 
be considered valid to the mouth of the Thames, where the amplitude of tide 
varies from about 5 ft at neaps to oft at springs, and the mean depth is of the 
order of 30ft. But a study of Southend residuals for the years 1928 to 1938 
eliminates even this possibility, for there are many instances on record when local 
winds were slight or even in the wrong direction for them to account for the 
observed interaction. 

It was therefore concluded that the interaction was generated when surge and 
tide were propagated from the open sea into the estuary under locally quiet 
meteorological conditions, and that attempts to compute the mode of propagation 


by numerical methods could possibly lead to a better understanding of the pheno- 
menon. 


3- Choice of numerical method 


To be satisfactory the method must be able to deal effectively with the principal 
non-linear terms in the equations of motion and continuity for one-dimensional flow 
along a channel. In addition a simple and direct control of the numerical input 
for any problem is highly desirable. Doodson (1956) has described his approach 
which consists essentially of integrating the fundamental equations, when trans- 


formed into finite difference equations, in distance along the channel. ‘This 
suffered from the disadvantage that the numerical input of any problem had to be 
in the form of a prescribed elevation at the head (barrier) of the estuary, and 
iterative methods were then necessary to find an appropriate prescription to pro- 
duce the required elevation at the mouth. This proved to be a reasonable procedure 
for the problem chosen, involving a depth of 128 ft, but to carry it out for a depth 
of 32ft was found to be impractical due to the generation of large shallow water 
tides. 

It so happened that at the time the choice of method was under consideration 
the author learned that numerical computations of a somewhat similar nature, but 
for an entirely different purpose, were in hand for an idealized mathematical model 
of the Thames. The problem was to determine the effects of a moveable surge 
barrier on river levels downstream of the barrier site, and this has now been 
successfully solved by Otter & Day (1960) using a modification of a method 
originally devised by Hansen (1956), and employing an electronic computer. 
The essential difference between the methods of Doodson and Hansen was 
that the latter’s consisted of numerical integration of the fundamental equations 
in time. 

The potentialities of this method are so great that it was decided to use it for 
the interaction problem; apart from differences in notation, the development given 
in the next section is substantially the same as that given by Otter and Day. The 
mathematical model they used was that of a channel of constant depth and expo- 
nentially varying breadth. This is a reasonable representation of the ‘Thames and, 
most important from the point of view of their preliminary theoretical investigations, 
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led to a considerable simplification of their formal solution for tides subject to 
linear bottom friction in the particular case of the Thames (see Appendix). In 
what follows an exponential channel of constant depth has also been used. 


4. The differential equations and their manipulation 
We shall denote by: 


h 
c 


W 


the mean depth of water, assumed constant; 

the elevation of the surface at time ¢ and position x, measured positive 
when above the mean level; 

h+g, the total depth of water; 

the time; 

the distance along the channel, measured from the mouth; 

the velocity of the water, assumed uniform from surface to bottom, and 
positive in the direction of increasing x; 

bWu = the quantity of water flowing past a section of channel per unit 
time at time f; 

bo exp(—2ax) = the breadth of the channel, where a is a constant; 

ne 
mr 
a coefficient of bottom friction; 

the acceleration due to gravity; 

the speed of a harmonic motion, i.e. 27/period. 


where m, n are integers; 


The usual differential equations may be written in the form: 


continuity : 


@0 at 
hey ad (1) 


motion: 


oP . (2) 


The non-linear inertia term u(0u/dx) has been omitted from the equation of 
motion since it complicates the computations out of all proportion to its importance 
in the particular model chosen. This was confirmed by preliminary calculations 
using a desk machine. 

If the estuary be partitioned by N+ 1 equidistant sections, commencing at the 
mouth (m = o) and ending at the barrier (n = N), and the origin of time be given 
by m = 0, the finite difference form of (1) at section m and time (m— 4)r is 


(Wm—Wm-s)n 1 Oma, n+a— Qm-4, n-4 
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by replacing Wm_;,n,, by the average of Wm,» and Wm-1,,1, and 
Wm-;,n-; by the average of Wm.» and Wiy-1,n-1. 

Then 


W’ n,n a4 Wm-i,n+Y¥m,nW m-1,n—1—8m,nWm-1,n+1 
where 


W’ = W1—y+8), ymn = dptm—y.n-1 Smn = $qUm-4,n+4- 


At section n+4 and time mr the finite difference form of (2) becomes 


uu 
= Wn Wan) “8 ! 


(Um+) —Um-3)n+4 ve 





m,n+} 


If Wmns4 is replaced by the average of Winn and Wmni1, and (u|u|)mns4 
bY Ums3,ns4|Um-4, n43| it follows that 


> 


, &t 
U' m+i,nty = Um-yntyt+— (Wan Wm,n+1) 
€ 
where 
ui = U1 +), &mignty = Wm,nt4lUm—j,ntal> 


6 
Wm,n+4 = 2kr/(Win,n+ Wm,n+1)- ( ) 


Inspection of equation (3) shows that if the values of W are known at all 
sections of the channel at time (m— 1)r, and the values of u are known at all points 
midway between sections at time (m—4)r, then W can be computed for all m for 
time mr. Similarly, equation (5) enables um, to be calculated for all points midway 
between sections solely from a knowledge of um—, and Wm. It is to obtain this 


state of affairs that the approximations leading up to equations (3) and (5) are re- 
quired. 


5- Boundary conditions 


At the barrier (n = N) the required condition is that Oy = 0 for all m. This 
affects the use of (3) at the barrier. 


In finite difference form the condition may be written 


(6uW)n+, = —(buW)y-, 


and hence (3) becomes 


W'mw = Wma,nt+2yWm-,n-1 (7) 
where 


W'y = Wy(1—2y), Ym,N = tPUm-;,N-1- (8) 


6. Arrangement of computations for desk machines 


To commence the calculations some arbitrary values must be given to W and 
u, and in the absence of an approximate solution to the problem to be considered 


the estuary water may be taken to be everywhere at rest with a depth correspond- 
ing to low water at the mouth, i.e. 


Um+, = 0 at m=o forall n+4 
Wm = W* at m=o forall n. 
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For all subsequent m, W* will be prescribed at n = 0 only, and represents a known 
disturbance of elevation ¢ travelling into the estuary from the sea. 

Two simple schedules are all that are necessary for the work to proceed. The 
first is in the form 


Ym 5m, W'm, Wm foreachof n=0 to n= N, 
and the second is in the form 
U'm+i, Om; Xm+y,Um+y foreachof n+}(m =o to N-1). 


With the initial arbitrary conditions listed above, schedule 1 is completed for 
m = I, using equations (3) and (4) for all sections except the barrier, for which 
equations (7) and (8) are used. Schedule 2 is then completed for m = 1 using 
equations (5) and (6). The calculations thus proceed alternately for W and wu until 
a solution is reached. For the problem of a continuous oscillatory input of constant 
amplitude and period the calculations must be continued over a sufficient number 
of cycles for the effect of the arbitrary initial conditions to have died away. In a 
typical example (Figure 3) only a little over 2 cycles of calculations were required. 


hours 
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Section 1-5 
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Fic. 3. Computed tidal elevations (a) and currents (b) at sections along 
an exponential estuary. ————— first cycle of computations, 
second cycle. Prescribed input at the mouth is 3 cos 30°F. 
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The calculations present no difficulties to experienced computers using desk 
machines, but suffer from the disadvantage that they are not amenable to summa- 
tion checks and need to be checked by repetition. The work is therefore most 
suited to an electronic computer, and it was programmed for the DEUCE machine 
belonging to the Computer Laboratory of the University of Liverpool by Mr 
J. M. Watt of that Laboratory, using a tabular interpretive programme. The 
author’s thanks are due to Mr Watt and to the Director of the Laboratory for their 
assistance in this aspect of the paper. 


7. Definition of the mathematical model; numerical data 


The average cross-sectional depth of the Thames between Richmond Lock 
and Southend ranges from 6ft at Richmond Lock to 4oft at Southend. A 
reasonable mean value of 


h = 32 ft. 
has been taken. 
Table 2 gives the mean breadth of the Thames at various points along its 


length, from which values of 2a in the expression b = bo exp( — 2ax) have been 
deduced. The value chosen for 2a in the model is 10-5 ft~} 


Table 2 


The exponential law of breadth for the Thames 
b = bo exp( — 2ax) 


Distance 
Place Ax = 


10°ft 
The Naze 
Southend 
Tilbury 


London 
Bridge 


The need to specify the elevation at the mouth of the estuary effectively results 
in the suppression, at the mouth, of waves reflected outwards from the estuary. 
In the case of a purely tidal problem the specification can be assumed to include 
any such reflected wave, and so far as the fundamental semidiurnal tidal oscillation 
is concerned there is no mis-matching problem. The quarter-diurnal and higher 
species of tide are generated in the estuary proper, however, and it is to be ex- 
pected that there will be some propagation of these species seawards. In the case 
of a transient oscillation such as a surge the mis-match can also be of importance, 
and indeed this difficulty is encountered in any method of injecting an input at the 
mouth of a tidal model, be it numerical, analogue or a small scale physical model. 

For the purposes of this investigation the mis-match has been minimized by 
taking the estuary twice as long as it is in nature. The “mouth” of the model 
estuary is thus somewhere in the vicinity of Lowestoft, the total length from mouth 
to barrier being approximately 117 miles. To a first approximation the speed of 
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propagation of a free wave can be taken as (gh), i.e. 32ft s~' Thus for a dis- 
turbance to travel from Southend to the barrier, and back to the mouth where it is 
reflected by mis-matching, and arrive again at Southend, a total distance of 234 
miles, would take some 11 hours. This interval of time is so long that mis-matching 
effects cannot materially influence the conclusions drawn from the computational 
results. A further safeguard in this respect is the decay experienced by any dis- 
turbance travelling seawards; an example of this effect is referred to in Section 
12 and illustrated in Table 6 for a semidiurnal tide of period 12 hours. 

Benefiting from the experience of Otter and Day, to ensure stability in the 
calculations the ratio «/r was kept above the greatest value of {g(h+)}* found in 
the channel, i.e. approximately 39fts~!. For 


tT = 600s, 2a = 10-5ft-!, N = 24 and k = 0-002 


this led to 
: &t 
€ = 0:2560x 10°ft, kr = 1-208, — = o-750s8"1, 
€ 
kp = 0°013320sft.-}, 
$q = o-or0311sft.}. 


Southend corresponds to section nm = 12 in the model, London Bridge to 
section 18 and Richmond Lock to section 24. 


8. Calculations for the tide alone 


Preliminary calculations by trial and error suggested an input at the mouth of 
3 cos 30°mr, which led to the results shown in Figure 3. These are in better 


agreement with observed neap tides in the Thames than might have been ex- 
pected from the crudeness of the model. The general shapes of the curves closely 
resemble those found in nature. Table 3 lists the observed ard model values of 
both range and mean tide level, and it will be seen that the largest differences occur 
near the barrier. The most likely explanation for these differences may be found 
in the fact that the bed of the Thames slopes upwards quite appreciably from the 
mouth, and Airy (1842, §343) showed the importance of such a slope in raising 
upstream mean tide levels above that of the sea. By taking the x-axis parallel to 
the bed of the estuary, and hence introducing a small component of the force of 
gravity into the equation of motion, it would seem probable that a closer approxi- 
mation to the tides in the Thames could be obtained in the model. For the pur- 


poses of the present investigation the amount of agreement obtained is all that 
could be desired. 


Table 3 
Comparison of observed and computed tides 


Mean tide level (ft) 
Range (ft) referred to that at Southend 


Section Computed Observed Computed Observed 
at neaps at neaps 
12 (Southend) 
18 (London 
Bridge) 
24 (Richmond 
Lock) 


o’7 


44 
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9. Comparison of computer results with those for the linear law of friction 


Doodson (1956) showed the usefulness of the linear law of friction if it is re- 
quired to study the propagation of a diurnal surge in the presence of a semidiurnal 
tide, when the shallow water distortions are not an object of interest. But his model 
was uniform in cross-section and rather deep (128 ft), and it is of considerable 
interest to learn whether the linear law is of practical value in much shallower 
water. In addition such an investigation, if successful, affords a welcome check 
upon the method of computation. 

The linear law expresses the frictional term in the equation of motion as fu. 
Comparing this with the quadratic expression ku|u|/h Proudman (1953, §151) 
showed that for a sinusoidal current 


f 37h . 
approximately, where U is the amplitude of the current. 

The theoretical solution for a semidiurnal tide in an exponential estuary of the 
dimensions used in this paper is given by equations (10, 11) in the appendix. 
Given the mean value of U along the channel all the necessary parameters, to- 
gether with the constant R, can be evaluated. From the computer results (‘Table 4) 
an average value for the estuary is 


U = 2:20 fts-1. 
Other numerical values required are 


c= 30° ht — 1°4544 x 1045-1, a= 05x 10-5ft-1, _ 32ft 


o2 
c2 = 1024 10-3 ft2s-2 Pm = 00434 x 10710ft-2 


fe = 0253 x 1075 0°328 x 10-5 
yr = —0o-172x 1075 0°828 x 107-5 
@ = 73°0 —34°2 
Values of amplitude (H) and phase lag (g) of the semidiurnal constituent of 
elevation in the DEUCE results were obtained by harmonic analysis and are given 
in Table 4, together with U and y for the semi-diurnal current. The origin of time 
is taken as high water at the mouth. Values of H’, g’, U’ and y’ derived directly 
from linear theory are also given in this Table. To compare them with the DEUCE 


results the phase lags have to be increased by 15°-6 to bring them to the same time 
origin, and the amplitudes have to be multiplied by the ratio 


R = (Ufrom DEUCE)/U’, or 9°87 x 10°. 


The last columns in Table 4 contain these adjusted values, and the agreement be- 
tween these and the DEUCE results is extremely good, proving the reliability 
of the linear law of friction for studying the propagation of the principal tidal 
constituent in a river resembling the ‘Thames. 


10. Calculations for the surge alone 


Whichever of the non-linear terms are responsible for the interaction it is plain 
that the greatest effect will be experienced when the total depth of water is least. 
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Negative surges were therefore the first to be investigated; if a negative surge at 
low water were to show no interaction effect, further researches along the pro- 
posed lines would be fruitless. 


The form of surge chosen was that of a single sinusoidal oscillation given by 


S = 2(cos15°t—1) for t = oto24, 
S=o0 for t < oandt > 24, 


where t is in hours. 


In a non-tidal estuary this surge is propagated as shown in Figure 4. The 
main points of interest are: 


(a) the speed of propagation of the first fall in level is in agreement with the 
theoretical celerity of 32 fts~}; 
(b) the small amplification as the surge progresses towards the barrier, 
and (c) the extremely rapid transmission of the lowest level. 
Calculations for this type of problem have been previously carried out by Reid 
(1957). 
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Fic. 4.—Computed surge elevations in a non-tidal exponential estuary. 
Prescribed input at the mouth is 2(cos 15°t—1) for t = 0 to 24. 


11. Calculations for tide plus surge 


The first combination of tide and surge considered was a 4ft negative surge, as 
defined in Section 10, phased so that its trough occurred at tidal low water at the 
mouth. Upon subtracting from the computed elevations the tide as computed in 
Section 8, the residues are found to differ appreciably (see Figure 5) from those 
of Section 10. In fact the computed surge exhibits precisely the interaction effect 
that is observed in nature in the Thames (see Figure 2). From Southend (Section 
12) to the head of the estuary the surge troughs occur on the rising tide, but the 
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development of this tendency can be noted as near the mouth as Section 6. More- 
over, it is remarkable that the trough of Figure 5 at » = 12 is 25 per cent lower 
than that of the surge alone (Figure 4), in accord with the observational result of 
Section 2. 
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Fic. 5. Computed negative surge elevations for case of input trough 
coincident with tidal low water. H and L denote tidal high and low 
waters respectively. 


To obtain complete confirmation that the model reproduces the observed in- 
teraction requires similar calculations to be made for other phase relations between 
tide and surge input, and the following cases have also been evaluated: 


tidal amplitude surge height peak at for residuals 
(ft) (ft) see Figure 


3 +4 high water 
low water 
rising tide 
falling tide 
as for +4 ft 

surge 
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Figures 6 and 7 exhibit the same phenomenon, viz., that in the estuary proper 
a surge is accentuated most on the rising tide. The computations for a 6ft surge 
have not been illustrated as they simply repeat, on a larger scale, the features 
mentioned. 

The diagrams also afford confirmation of one of Proudman’s conclusions 


(1955b), that “.. . the surge up the estuary is less at the time of high water than 
it is at the time of low water ...’’. 


12. A simple physical explanation for interaction in the Thames 


In Figure 8 are drawn the computed elevations of tide plus surge in the case of 
a negative surge coincident with tidal low water at the mouth; this is the case for 
which the residues are shown in Figure 5. In the same figure are shown the re- 
sults of arithmetical superposition of the computed surge in the non-tidal estuary 
(see Section 10) upon the computed undisturbed tide of Section 8. If we consider 
the difference between these two sets of curves to be due to the influence exerted 
by the surge upon the propagation of the tide along the estuary it is apparent that, 
with respect to the undisturbed tide, the disturbed tide becomes increasingly later. 
We may therefore conclude that the presence of a negative surge causes the re- 
tardation of the tide up the estuary. The converse to this, that a positive surge 
hastens the tide, is supported by the contents of Table 5, which lists the time 
shifts of maximum level for 8 combinations of tide and positive surge. ‘Two time 
shifts are given for each combination, as the figures are derived from computations 
covering one diurnal surge. 

A fall in level, i.e. a negative surge, can retard the tide in two ways. The rate of 
progression of a free wave will be reduced by a reduction of depth; a positive surge 
will have the opposite effect. Bottom friction also has the effect of retarding a 
wave, and as bottom friction is inversely proportional to depth, a negative surge, 
by reducing the depth, will increase the friction and retard the wave even further; 
a positive surge will have the opposite effect. 

Thus there are two independent mechanisms whereby a surge can influence 
the time of arrival of a tide, and they provide an extremely simple qualitative 
explanation for the preference both positive and negative surges show for the rising 
tide. Figure 9 shows a simple tide curve (A) displaced by half an hour (B). Let 
A represent the tide under normal conditions, and B the retarded tide due to a 
negative surge. It will be seen that the only state of tide which will contribute a 
further lowering of level is the rising tide. Now consider the case of a positive 
surge hastening the tide; this time B will represent the normal tide and A the 
hastened tide. Once again it will be seen that the only state of tide which will 
contribute a further appreciable rise in level is the rising tide. 

One possible objection to the above explanation might be that the tides in the 
Thames are more nearly in the form of a standing oscillation than a progressive 
wave. This opinion has been very largely accepted in the past on the assumption 
that the estuary is uniform in cross-section. Otter & Day (1960) and Abbott 
(1960) have now shown that an exponential cross-section favours the progressive 
wave. 

The theoretical solution for the principal semidiurnal tidal constituent, given 
in the appendix to this paper, clearly indicates two progressive waves for the chosen 
exponential channel. The incoming wave is given by the first expression in each 
of equations (10) and (11), and the reflected wave by the second expression. For 
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the incoming wave it is shown in Section g that the current leads the elevation by 
6(= 73°) and that the amplitude of elevation increases from the mouth as exp(rX) 
with r negative and X = o at the barrier; in the reflected wave the current lags 
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Fic. 6.—Computed positive surge elevations for cases of input peaks 
coincident with tidal high water ( ) and tidal low water 
(—__———). H and L denote tidal high and low waters respectively. 


on the elevation by 7—¢(= 214°) and the amplitude of elevation decreases 
towards the mouth as exp(—sX) with s positive. Table 6 gives values of exp(rX) 
andfexp(— sX) at points along the channel, and it is evident that the reflected 


wave is insignificant at the mouth, supporting the views of Otter & Day and of 
Abbott. 
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Fic. 7.—Computed positive surge elevations for cases of input peaks 
coincident with the rising tide (—————) and the falling tide (- - - - - -). 
H and L denote tidal high and low waters respectively. 
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Table 5 


Time shifts of maximum level relative to the tide alone, for various combinations of 
tide and positive surge. Amplitude of tide is 3 ft, and units of time are minutes. 


Surge peak at mouth coincident with 
High water Falling tide Low water Rising tide 
Lags on tide alone 
Surge of 4 ft 
3° —35 
15 ? 
14 —4I1 


7 —44 
— 1 —46 


Surge of 6 ft 
° : 35 —50 60 —55 
—10 25 —50 40 —65 
—21 14 —51 34 —66 
—39 —- —59 26 —89 
—4!I — I —66 24 —96 


Lags on tide alone relative to those at Section o. 


Surge of 4 ft 
—I5 ? —10 
—16 — 6 —16 
“Mo + 8 
—31 —11 —21 


Surge of 6 ft 


—10 ° —20 
—21 — I —26 
ae ee —34 
—36 —16 —36 


Table 6 


Relative amplitudes of incoming (e"*) and reflected (e~**) waves for an exponential 
estuary resembling the Thames 


r= —0o'172 1075, s = 0828 x10°5 


Section erx e~tx 
o (Mouth) 0°347 0006 
6 0°453 0°022 
12 0°59°0 0°°79 
18 0-768 0280 
24 (X = 0) 1000 1-000 
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Accepting the progressive wave regime of the tides in the Thames in personal 
correspondence, Proudman points out that formula (42) of his 1955a paper may 
be used to show that: 

“For equal sequences of meteorological conditions, the maxima of apparent 


surges are greater when they occur at a quarter tidal period before tidal high 
water than when they occur at a quarter tidal period after tidal high water”. 


This statement is in accord with the results of the present paper. 

In passing, it is worth noting that Doodson’s arguments of 1928 (see Section 
2) are not so very different from those propounded here; he suggested a modi- 
fication of tidal range by a surge whereas I suggest a time shift; he suggested wind 
stress as a potent factor whereas I suggest bottom friction. 


hours 
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iz 
Fic. 8.—Comparison between computed elevations for an input of tide 


plus negative surge (—————-) and the superposition of computed 
elevations for tide alone and surge alone (— - - - - - ). 
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The above arguments relate to the time shift of a tide in the presence of a surge. 
Similar qualitative conclusions may be arrived at if the tide is considered as 
shifting the surge in time relative to the surge in non-tidal waters. To attempt a 
determination of the relative importance of these two aspects of the problem would 
appear to be profitless. What clearly emerges from the foregoing is that “‘surges’’ in 
the Thames, as defined by subtracting the undisturbed tide from the observed 
levels, cannot be forecast solely in terms of meteorological variables. Indeed, fore- 
casting formulae for the Thames should be ideally designed to predict, not a surge 
which can be added to the predicted tide, but tide plus surge as a whole. 














Fic. 9.—A semidiurnal oscillation (A) displaced half an hour (B). 


The computer results, when viewed in the light of the physical mechanisms 
suggested, indicate that interaction is most likely when tide and surge have been 
travelling together for some distance in quite shallow water. This would explain 
why no interaction has been found at Tynemouth, Immingham and Lowestoft, 
the mean depth of water being too great. By the same token it would be expected 
that interaction should be found well up the River Humber, in the Wash and in the 
River Severn. 


13. The propagation of maximum levels 


It was stated in Section 1 that some practical results were required from the 
present investigation; in particular there is a need for a simple method of includ- 
ing the effect of interaction on the propagation of dangerously high levels into 
and along the Thames. 

Table 7 lists the maximum levels reached at sections 0, 6, 12, 18 and 24 as 
obtained from the results of 8 combinations of positive surges with the same 
tide. Two heights are given for each combination as the figures are derived from 
computations covering one diurnal surge; the heights are referred to the undis- 
turbed level of 32 ft. In the same table are given the amplifications of these levels, 
as they travel up the estuary, relative to the level at the mouth. They show that 
the amplification is a function of the level at the mouth, the nature of the function 
being illustrated in Figure 10. It is clear that there is a strong tendency for the 
higher, more dangerous levels to be amplified less than the lower levels. Another, 
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but less comprehensive way of putting this, is that dangerous “‘surge’”’ peaks tend 
to avoid tidal high water. 

Whilst of considerable practical importance, the amplifications given in Figure 
10 cannot be applied directly to the ‘Thames, as the mathematical model cannot 
be considered a sufficiently close approximation to nature; the results nevertheless 
indicate profitable lines for future research into both mathematical models and 
observational data. 


Table 7 


Maximuin levels of tide plus surge referred to the undisturbed level of 32 ft. Amplitude 
of tide is 3 ft, and units of height are ft 
Surge peak at mouth coincident with: 
High water _‘ Falling tide Low water Rising tide 
Maximum levels 

Surge of 4 ft Tide alone 

Section o . ° 3°7 5°2 : : 6°5 3°0 

6 e . 4°6 6°3 . . 76 4'1 

12 6-0 8-1 ‘ , 9°3 5°7 

18 . 8-4 10°8 ‘ 11°6 8-2 

24 : 9°9 12°7 130 10° 


Surge of 6 ft 
30 8-2 4°1 64 6°4 4'1 8-2 
4:0‘ 96 5:0 7.39 52 9°5 
5‘2 115) 663, 94 88 6°7 
747 144 85 I2"l I1‘o 90 
9°4 16°5 10°0 14'2 123 10°3 


Ratios of maximum levels to those at Section o 
Surge of 4 ft 
1°33 1°20 1°24 I'2r 1'Ig 1°30 
180 1°48 1°62 1°56 1°46 1°70 
2°63 189 2°17 2°08 1°94 2°35 
3°20 2°18 2-68 2°44 2°58 2°73 


Surge of 6 ft 
1°33 1°22 I'lg 1°17 1°27 
1-71 1°54 1°47 1°38 1°63 
2°57 2°07 1°89 «61°72 2°19 
3°13 2°46 2°22 «1°92 2°51 


APPENDIX 


Propagation of tides in a channel of constant depth and exponentially 
varying breadth under the influence of the linear law of friction. 


We shall use the same notation as in Section 4 except that it is more con- 
venient to take X = L—x where L is the length of the channel, so that X = 0 at 
the barrier. Hence b = bo exp(2aX), and the positive direction of u is towards the 
sea. Also let bu = V. The equations of continuity and motion for first order 
terms only are 


hVx = —bl,, bglxy = —Vi-fV 
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where suffixes X, ¢ represent partial differentiation with respect to these variables. 
For a single harmonic constituent of speed o, V is given by the real part of 
exp(iot)F(x), and elimination of { from the above equations gives 

of 
Vxx—2aVx¥+— (: —1 -)V = 0. 
c o 
The solutions of this homogeneous differential equation are derived from the roots 
of the auxiliary equation 


pt—zap+ = (1-17) =O (0) 


and are given below for the condition of no current at the barrier (u = 0 at X = 0). 


Maximum teveil ef mouth (ft) 
6 Zz b , 
T t tT ia 
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Fic. 10.—Amplification of maximum levels up an exponential estuary 


relative to those at the mouth, as a function of the maximum level at the 
mouth. 


O for tide+surge; @ for tide alone. 


Case I. f = 0. Auxiliary equation: p* — 2ap + o7/c? = o 
(1) a > a/c. The roots « and £ are real and unequal, 


a= at(at—o8/c)), B = a—(a*— oct) 
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and the solution is 


{ = Ris exp(rX)—rexp(sx)] cosot, u= — , R{exp(rX)—exp(sX)] sin ot 


where 
r=a-—2a, $s = B—2a,and Ris an arbitrary constant. 
(2) a = o/c. The roots are real and equal, and the solution is 
C = Rexp(—aX) (1+a@X) cos ot 


“= ~ Rexp(—aX) X sin ot. 


(3) a < o/c. The roots are complex, « = a+ir, 8 = a—ir where, r? = o2/c2—a?, 
and the solution is 


{ = Rexp(—aX) (rcosrX+asinrX)cosot, u= > R exp(—a@X)sinrX sin ot. 


The particular case of a uniform frictionless channel arises here with a = 0, 
giving the well-known solution 


" Peee ft te oO 
¢ = R-cos—X cosot, “= R>sin— Xin ot. 
” ape c 


Case II. f # 0. 
(1) a > o/c. The roots «, B of the general auxiliary equation (g) are complex: 


a,8 = at(«+ip) 
where 


K2—py2 = a2?—o/c2?, = 2xp = affc?. 
The solution is 
C = R{(u? +52)! exp(rX) cos(ot + wX — 6) + (u? +1?) exp( —sX) cos(ot —uX —¢)} 


(10) 
“— -; R{exp(rX) cos(ot + .X) — exp(—sX) cos(ot — uX)} (11) 


where 
r= Kk-—4@, s= «+a, 6 = tan“! s/p, @ = tan“! r/p. 


(2) a < a/c. The roots are given by «, 8B = a + (u+ix) and the solution is 
obtained by interchanging « and u in (10) and (11). 

The particular case of a uniform channel arises here with a =o and 
p2— x? = —o*/c?, for which the solution is 


C = R(p?+ x2)! {exp(X) cos(ot + xX — 0) + exp(— 4X) cos(ot — «X — 6)} 
u=— R {exp(uX) cos(ot + «X) — exp( —X) cos(ot — «X)} 


with 
6 = tan—lu/«. 
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(3) a = o/c. Here « = x = (of/2c®)*, and the solution is the same as either of 
cases II(1) or (2). 


Remarks on solutions 


In the absence of friction and geostrophic effects (Case I) the incoming and 
reflected waves for all species of tide combine to give the well-known equations 
for standing oscillations. The introduction of friction (Case IT) affects the waves 
so that in each of them the elevation and the current are not in phase, and standing 
oscillations cannot result from their superposition. In case II(1) the speed of 
propagation of both waves is o/u, in case II(2) it is o/«. With the dimensions 
chosen for the model, 


ac = 16x 1074s"! 


and it follows that Case II(z) is applicable for all tidal constituents of period greater 
than 2z/ac, or 10-9 hrs. 

One particular case of interest arises under II(1) when x = a. Then yu = a/c, 
f = 2ac and 


{ = R{(u2 + 4a?) cos(ot + .X — 0) +p exp( — 2aX) cos(ot — pX)- 


= -; R{cos(ot + 1X) —exp(— 2aX) cos(ot —2X)}. 


It will be seen that the incident wave is of constant amplitude, both in elevation 
and current, and is propagated with speed c up the estuary. The reflected wave 
decays at the rate e~2¢* as it travels seawards, and could be negligible at the mouth. 
This is the special case used by Otter & Day in their analytical work on the Thames 
barrage, and first discussed by Harris (1907, 5, p. 295). 
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Heat Flow in the Austrian Alps 


Sydney P. Clark, Jr. 
(Received 1961 March 21) 


Summary 


Data on underground temperature obtained during the construction 
of the Arlberg and Tauern tunnels in Austria have been combined with 
measurements of the thermal conductivity of 42 samples of rock from 
near the tunnels to calculate the terrestrial heat flow. The value in the 
Arlberg is found to be (1-9 + 0-2) x 10~® cal/cm? s, and in the Tauern, 
(1-8+0-2) x 10-* cal/cm? s. The new results are in good agreement 
with the value 1-9 x 10-6 cal/em?s found earlier in the Loetschberg 
tunnel in Switzerland, and indicate that relatively high geothermal 
fluxes extend into the eastern Alps. The high flux can be attributed 
to radioactive heat generation in a thickened crust. 


1. Introduction 


Determinations of terrestrial heat flow in the Alps, based on the temperatures 
measured in the three longest Swiss railroad tunnels, were reported in a previous 
paper (Clark & Niblett 1956; hereafter referred to as 1). In the summer of 1960 
the writer had the opportunity to visit two of the great Austrian tunnels and collect 
samples for the measurement of thermal conductivity. The results, coupled with 
the temperatures observed during construction of the Arlberg and Tauern tunnels, 
yield two determinations of heat flow in the eastern Alps, which provide an 
interesting extension of the Swiss data. 

The locations of the Arlberg and Tauern tunnels are shown in Figure 1. 
The Arlberg runs through gneisses and schists of the pre-Mesozoic basement 
beneath the Arlberg pass, which connects the Tyrol with the Vorarlberg. The 
tunnel was completed in 1884. The Tauern tunnel runs through granites and 
granite gneisses which out-crop in the so-called Hohe Tauern Window. It lies 
about 220km to the east of the Arlberg, and it was completed in 1908. 


2. Temperature measurements 


Temperature profiles are shown in Figures 2 and 3. In the Arlberg the tem- 
peratures of the rock were measured with mercury thermometers in bore-holes 
o°8 m deep in the walls of the tunnel. The measurements are described by Wagner 
(1884), from whose paper the data were taken. The temperature of the air of 
the tunnel was measured simultaneously; it was usually greater than the rock 
temperature, but the difference rarely exceeded 3°C. Judging by the irregularities 
in the temperature profile of Figure 2, the measurements are uncertain by about 
1°C. 
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The temperatures in the Tauern tunnel were taken from handwritten reports 
of the construction company, which are preserved in the Verkehrsarchiv in Vienna. 
The technique and accuracy of the measurements are presumably roughly the 
same as in the Arlberg. The temperature measured at the point 4km from the 
north portal appears to be definitely low in relation to the surrounding measu re- 
ments, however, and a value interpolated from adjacent observations was used for 
this point (dashed line, Figure 3). 

The Tauern tunnel encountered springs between 6 and 6-5 km from the north 
portal. These appear to have cooled the rock in their vicinity, and only the tem- 


peratures measured within 5-5km of the north end of the tunnel were used to 
determine the heat flow. 
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Fic. 1.—Location map. The Arlberg, Tauern, Gotthard, Loetschberg, 

and Simplon tunnels are shown as heavy dashes and are identified by 

the first letters of their names. The numbers following the letters are the 
heat flows in units of 10~* cal/cm?s. 


The temperatures at the portals of the tunnels were estimated by extrapolation 
of the temperatures within the tunnels. In the Arlberg a value of 6-3°C was 
obtained in this way, which is reasonable for the elevation of the tunnel. A more 
complicated situation is encountered in the Tauern, since near the portal the 
tunnel runs in valley fill with only a few tens of metres of cover. About o-5 km 
from the portal it passes under the steep cliffs which form the south wall of the 
Anlauftal; the temperature at this point, found by extrapolation to be 4-7°C, 
is taken to be the same as the temperature at the portal. This temperature is 
unusually low for the elevation of the tunnel, probably because the point of 
observation lies at the foot of a steep cliff which faces northward. 

Neither tunnel is level, and it is convenient to correct the temperatures in 
each to a common elevation. This was done using the uncorrected gradients 
given in Table 2, which were determined by expressing the observed temperatures 
as functions of vertical cover without correction for topographic irregularities. 
The Arlberg was corrected to an elevation of 1 300m; the maximum correction 
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Fic. 2.—Topographic section and temperature profile, Arlberg tunnel, 
Temperature : ————-; elevation: 
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Fic. 3.—Topographic section and temperature profile, Tauern tunnel. 
Temperature : ————; elevation: — ---— -. 
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was 0:2°C. ‘Temperatures in the Tauern were corrected by 0-5°C at most to 
bring them to an elevation of 1 200m. 


3. Geothermal gradients 


Topographical corrections were made for stations at intervals of o-5_km from 
the portals of both tunnels. The correction was the same as used in I (see also 
Clark 1957). Account was taken of the evolution of the topography by assuming 
that an initially plane surface was uplifted by an amount L in time ¢ and simul- 
taneously eroded to the present topography. The present (projected) position of the 
initial surface is a distance D above the tunnel, and the depth of a station is denoted 
by z. Calculations were made for t = 25 x 108 years, L = 4km, and D = 3-5km, 
and for t = 4x 10% years, L = 3km, and D = 2:4km. These values were also 
used in I, and hence a direct comparison of the results is possible. 

Corrected results are given in Table 1. Symbols not previously defined are «, 
the geothermal gradient; «’, the rate of decrease of surface temperature with 
elevation; and a, which stands for 2/(x«t)', where « is the thermal diffusivity. 
The topographic correction is contained in the term A(t), which is found from 
topographic maps in a manner similar to the correction in gravity surveys. 

From the measurements of thermal conductivity and density given below, « 
is found to be roughly 50km?/10® years. The gradient «’ was taken to be 4°67 
deg C/km, a value that was obtained from a compilation of soil temperatures 
briefly described in I. For comparison, Kerner (1905) found a gradient of 
4°4degC/km in the central Tyrol and elevation-dependent gradients of 5-5 
deg C/km at 1250m and 4-odegC/km at 2000m in the North Tyrol. His results 
were based on the decrease in temperature of springs with elevation. 

A plot of corrected temperature against corrected depth should yield a straight 
line with slope «—«’ and intercept 7», where To is the surface temperature at 
tunnel level. Parameters of least-squares lines, fitted to the data for the various 
assumptions about the physiographic history, are given in Table 2 and two cases 
are shown in Figures 4 and 5. In the Arlberg tunnel the points representing the 
stations near the portal fall systematically above the line, and the temperature at 
the portal falls about 1°C below it. No point departs from the least-squares line 
by more than 1-1°C, however, and hence the results are within the uncertainties 
of the topographic correction and the measurements of temperature. 

The temperatures at stations 1-0 and 1-5 in the Tauern tunnel fall as much as 
3°6°C below the least-squares lines, and the intercepts of the lines are 1-5°C or 
more beneath the temperature at the portal. Three plausible explanations of the 
low temperatures can be given. ‘The stations are close to the steep cliffs forming 
the south wall of the Anlauftal, and the topographic correction, which is approxi- 
mate, may fail in cases of such rough topography. The cliffs face to the north, 
which may cause the temperatures at the surface, and consequently the tempera- 
tures at depth, to be abnormally low. Finally, the rocks in this vicinity are thermally 
anisotropic, conducting best in the plane of their foliation. The strike of the 
foliation is roughly perpendicular to the trend of the valley, so that there is a 
tendency for heat to be channelled towards the cliffs. This could also lead to low 
underground temperatures. 

The data given in parts B and C of ‘Table 2 were calculated by retaining or 
rejecting stations 1-0 and 1-5. The difference in gradient is 1-5 deg C/km at most, 
or about 7 per cent. Without the low points, the slope of the curve depends 
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critically on the temperature at the portal since the remaining stations range over a 
corrected depth of only about 300m and the shallowest of them is 750m in depth. 
If the standard error of the slope is calculated in the usual way, no account is 
taken of the uncertainty of about 1°C in the temperature at the portal. The 
standard deviations of the slopes given in part C of Table 2 were found by deter- 
mining the variance of the slope as usual, and increasing this value by 1 (deg C/km)*. 
This is very nearly equivalent to assuming that the temperature at the portal has a 
standard deviation of 1°C which is independent of the scatter of the remaining 
points. The standard errors calculated in this way are close to those given in part 
B of Table 2. The differences between the g-point and the 11-point solutions are 
much less than the differences between the extreme assumptions about the topo- 
graphic history, and no serious additional uncertainty in the final results is intro- 
duced by the lack of an objective way of deciding which is superior. 


Table 2 


Parameters of least-squares lines 


t D ao" Intercept Standard 
10° years m deg C/km °C deviation 


"C 


A. Arlberg Tunnel 
oe) aa an 15°01 +0°94 7°24 o'72 
25 4000 3 500 13°10 +0°74 7°26 0°62 
4 3 000 2 400 12°03 t0°75 7°33 0°66 
Uncorrected 14°'O2 + 1°14 8-04 


o'9! 


B. Tauern Tunnel, including stations 1-0 and 1°5 
oa) len sb 21°75 +1°88 2°99 
25 4 000 3 500 18°20 +1°47 3°14 
4 3 000 2 400 16°47 + 1°42 3°21 
Uncorrected 13°97 1°39 5°23 


C. Tauern Tunnel, excluding stations 1-0 and 1°5 

ees one 20°20+1°44 4°80 
4 000 3 500 16°94 + 1°40 4°77 
3 000 2 400 15°31 +1°37 4°81 


4- Thermal conductivity 


Measurements of thermal conductivity were made on samples collected at 
the surface above the tunnels. Conductivity was determined in a modified divided 
bar apparatus, of the kind described by Birch (1950), at a mean temperature of 
25°C. The samples could be loaded with axial pressures up to 2 500lb/in?; the 
effect of this load was less than the experimental error, which is estimated to be 
about 2 per cent. 

The rocks of the Arlberg presented a problem in sampling, because they are 
steeply dipping garnetiferous phyllites and micaceous gneisses. It was difficult 
to collect specimens from which disks, 1-5 inches in diameter, could be cut with 
the plane of the disk perpendicular to the plane of the foliation of the rock. Eight 
samples were collected from outcrops along the road which winds over the Arlberg 
pass. Despite efforts to make the sampling as representative as possible, it is un- 
doubtedly biased in favour of the more resistant, less fissile units. ‘These are 
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Fic. 4.—Corrected temperatures as function of corrected depth, Arlberg tunnel. 


The stations are numbered by their distances in kilometres from the 
east portal. t = 25 x 10° years, L = 4km, D = 3°5km. 
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Fic. 5.—Corrected temperatures as function of corrected depth, Tauern 
tunnel. The stations are numbered by their distances in kilometres from 
the north portal. ¢ = 25 x 10% years, L = 4km, D = 3°5km. 
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probably the best thermal conductors, but since the structure is such that they are 
arranged in parallel with the poorer conductors, these beds are the most important 
to sample. 

Measurements of the conductivity in a direction parallel to the plane of the 
foliation were made on seven samples, and the conductivity normal to the foliation 
was measured for all eight specimens. Results are summarized in Table 3. Because 
of the steep dip of the foliation, the “‘parallel’’ specimens should be used to cal- 
culate the heat flow. These form a fairly homogeneous group with respect to con- 
conductivity, despite considerable differences in appearance and mineralogical 
constitution. Lacking a far more elaborate sampling procedure, we can do no 
better than calculate the heat flow from the mean of the seven measurements. 


Table 3 
Thermal conductivity 


No. of Thermal conductivity, 10~* cal/ems deg C 
measurements Orientation mean range 


Arlberg Tunnel, gneiss and schist 


parallel 10°74 +0°37 9°8-12'5 
perpendicular 7°23 +0°36 5°0-8'0 


Tauern Tunnel, granite and granitic gneiss 
Tauern 
13 parallel 8-87 +0°18 
9 perpendicular 6°85 to'10 
5 isotropic 7°50+0°02 


Mean density, Arlberg, 7 specimens, 2°686 +0°033. 
Mean density, Tauern, 19 specimens, 2610 + 0-003. 


The rocks in the vicinity of the Tauern tunnel are rather uniform granites and 
granitic gneisses. Samples were collected from fresh talus slopes along the Min- 
dener Weg, a footpath which crosses the mountains in a direction parallel to the 
tunnel and about 2km to the east. 

Most of the specimens showed clear-cut thermal anisotropy (Table 3). Accord- 
ing to Exner’s tectonic map of the region (Exner 1957, Plate VIII), the dip of the 
foliation ranges from slightly less than 30° to less than 60° in the vicinity of the 
tunnel. If the mean dip is 40°+ 10°, the vertical component of the conductivity 
of the anisotropic rocks is (7°68 + 0°37) x 10-%cal/emsdegC; most of the un- 
certainty in this figure arises from lack of quantitative information about the 
geologic structure between the tunnel and the surface. This number does not 
differ significantly from the mean conductivity of the few isotropic specimens. 

The mean temperature of the rock between both tunnels and the surface is 
about 10°C to 12°C. The observed conductivity should be increased by about 
2 per cent to take account of the difference between this temperature and the 
temperature of the measurements. We adopt the round figure (11-00 + 0°37) x 
10-8 cal/emsdegC for the conductivity of the rocks in the Arlberg tunnel, and 
(7:83 + 0°37) x 10-3 cal/ems deg C for the conductivity in the Tauern. 
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5. Flow of heat 

The flow of heat for the various assumptions about the topographic" evolution 
are given in Table 4. Two versions of the results for the Tauern are included; 
they correspond to retention or rejection of stations 1-o and 1-5. The results for 
the two tunnels are in good agreement with each other. At worst the difference 


only slightly exceeds the combined standard errors for the same assumption about 
the topographic history. 


Table 4 
Flow of heat 


Heat flow, 10~* cal/cm?s 


Arlberg Tauern Tauern 
II points 9 points 


co 2°'16+0°12 2°07 +0'18 1°95 tO°1S 
25 1°95 toO'10 1°79 +0'14 1°69 +0°14 
4 1°84 +0°10 1°66 +0°14 1'°56+0°13 
Adopted value 1'9 +0°2 1'8+0°2 


The Austrian values are also in good agreement with results for the Loetschberg 
tunnel given in I. They are rather lower than the value 2-2 x 10-6 cal/cm?s adopted 
in the Simplon and considerably higher than 1-4 x 10~® (corrected to 1-6 x 10-6) 
cal/cm?s, as found in the Gotthard tunnel. Estimates of the heat flow in the two 
Austrian tunnels were hazarded in I. They are in fair agreement with the present 
determination in the Arlberg, but in the case of the Tauern the estimate was too 
low. (t was thought that the estimated fluxes were, if anything, too high, but this 
does not prove to be the case. 

An estimate was also made in I of the heat flow in the Albula tunnel in eastern 
Switzerland. It was later discovered that the gradient in this tunnel was based 
on a single measurement of temperature made several months after construction 
was completed. Hence this estimate is wholly unreliable. 


6. Discussion 


On the basis of the estimates of the heat flow in Austria, it was concluded in I 
that the heat flow in much of the Alpine region was about 1-5 x 10-6 cal/cm?s. 
The higher values in the Simplon and Loetschberg tunnels were considered to 
represent a more or less local anomaly. The previous discussion requires modifica- 
tion in the light of the present results, which show that high values of heat flow 
persist far to the east of the earlier determinations. 

As in I, we shall assume a layered crust for purposes of discussion. ‘The crust 
outside the mountains is assumed to consist of a layer of granitic rocks, 20km 
thick, lying on 10km of basaltic rock. In the mountains the granitic layer is 
thickened to 35 km and the thickness of the basaltic layer is unchanged. If 5 x 10-18 
cal/em’s is produced by radioactivity in the granitic layer and 10~!%cal/cm$s 
in the basalt, the heat flow in equilibrium with crustal radioactivity is 1-85 x 10-6 
cal/cm?s in the mountains and 1-1 x 10~®cal/cm?s elsewhere. This model leads 
to plausible heat flow in the lowlands and to nearly sufficient heat flow in the 
mountains; the remainder can be assumed to originate beneath the crust. 
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If the root of the Alps was formed relatively recently, the heat flow at the 
surface may not be in equilibrium with the radioactivity of the root, and account 
of this effect may be taken by the approximation used by Birch (1950, p. 618). 
If the root was formed 25 x 10° years ago, it is found by this method that the 
heat production in the granitic layer must be increased to 6x 10~!8cal/cm*s in 
order to obtain a heat flow of 1-8 x 10-6 cal/cm?s from crustal radioactivity. The 
equilibrium heat flow outside the mountains is then 1-3 x 10~6cal/cm?s. Both 
this figure and the value of the radioactivity of the granitic layer are acceptable in 
view of present knowledge of these quantities. ‘There appears to be no difficulty in 
accounting for the high heat flow in the Alps on the basis of radioactive decay in a 
thickened crust. 

The relatively low heat flow in the Gotthard tunnel can be explained in an 
ad hoc fashion only. Heat production in the underlying crust may be unusually 
low in its vicinity. Values of 3 x 107!8 to 4 x 10718 cal/cm*s would be demanded, 
with the exact figure depending on the size and shape of the region of low heat 
production. Alternatively the vicinity of the Gotthard may have undergone less 
uplift in recent times than other parts of the Alps, despite the fact that this tunnel 
passes under the highest-standing part of Switzerland. At present there seems to 
be no way to decide between these explanations. 
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The Deviation of the Vertical at Bidston in response 
to the Attraction of Ocean Tides 


G. W. Lennon 
(Received 1961 March 24) 


Summary 


The periodic deviation of the vertical in sympathy with the flux of 
tidal waters is of particular interest in the study of Earth tides and 
associated phenomena. In order to assist Earth tide studies at Bidston, 
the attraction exerted by the presence of tidal waters has been evaluated 
for this station. The major lunar semidiurnal constituent of the 
phenomenon, (M2), has been obtained by direct computation from a 
knowledge of the tides in the vicinity, and a method has been devised 
for the inference of the principal solar semi-diurnal constituent, (So), 
and the larger lunar elliptic semi-diurnal constituent, (Ne), based upon 
regional values of the relationships between the phases and amplitudes 
of the appropriate constituents in the ocean tide. 


The use of the horizontal pendulum for the measurement of variations in the 
tilt of the Earth’s surface is particularly susceptible to the attractive effect of 
oceanic tides since deviations of the vertical occasioned thereby are recorded by 
the instrument as apparent tilts of the surface and consequently present confusing 
elements in the study of real crustal movement. The phenomenon is especially 
significant at Bidston which is sited on a peninsula between two tidal estuaries and 
only two miles distant from the shores of the Irish Sea in which are experienced 
tides up to 30 ft in range. Of the total recorded periodic tilt at Bidston, less than 
one tenth part can be attributed to the true earth tide and the remainder is due 
mainly to the effects of oceanic loading and oceanic attraction. In earlier work on 
crustal tilts at Bidston, the effects of loading and attraction have not been separated 
but have been considered as a unit having the same phase and azimuth. This 
assumption may well be justified within the limits of accuracy of the investigation, 
but it is felt that as the results of more detailed observation become available, 
allowing greater precision, then this assumption cannot be accepted without 
examination. It is possible for example that the loading effect may show a lag upon 
the loading tide whereas the attraction effect must necessarily be simultaneous 
with the tide. Furthermore, recent work suggests the fundamental influence of 
regional geological structure in determining the azimuth of crustal tilt so that, in 
this respect also, loading may be affected independently of attraction. Again, the 
range of influence of the two effects may be quite dissimilar in character thereby 
giving differences in both azimuth and phase. It is clear that in an attempt to 
understand the complex mechanism of periodic surface movement the evaluation 
of the part played by tidal attraction is a significant and essential step. 


64 
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1. The source of tidal information 


The primary requirement of an investigation of this type is a store of precise 
knowledge of ocean tides, which is not generally available, and it is fortunate that 
the waters adjacent to the British Isles have been the subject of careful investigation 
through the years, particularly at the Tidal Institute. As a result, the disposition 
of the principal lunar semidiurnal constituent, (Mg), in this area is known just 
sufficiently well for a direct computation of the Me constituent of tidal attraction 
to be attempted. In this respect the basis is the British Admiralty Co-tidal Chart 
of the Me constituent, which chart is itself founded upon work by Proudman and 
Doodson (1924) for the North Sea and by Doodson & Corkan (1932) for the 
English Channel and the Irish Sea. 

Nevertheless the scope of the existing work was not immediately adequate for 
the study in hand and it was necessary first to extend the co-tidal chart westwards 
so as to include the eastern part of the Atlantic Ocean and also northwards to 
Iceland and the Faroes. In this task the coastal directions of the co-tidal and co- 
range lines were first computed, wherever possible, according to the theory de- 
veloped by Proudman & Doodson (1924) and using the phase lags of the 
maximum tide and of the maximum flood stream. With the assistance of the 
azimuths so obtained it was then possible to fit co-tidal and co-range lines to the 
published tidal data derived from coastal observation at specific locations. Figure 1 

















Fic. 1.—Co-tidal chart for the Me constituent. 
co-tidal line showing the phase lag of M2 


— O—— co-range line showing the amplitude of Mz in feet (Roman 
numerals) 
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illustrates the general framework of the resulting chart although it should be 
noted that the tidal basis actually utilized in the calculations was defined in con- 
siderable detail and with a fair degree of precision. 


2. The computational method for the M; constituent of attraction 


The problem of the deviation of the vertical has been the concern of geodesists 
for a considerable period since the effect of local geology and relief upon the tech- 
niques of precise surveying has had to be considered. As a result methods of 
determination of the magnitude of the effect from the physical phenomena are 
well known and the selection of a technique is mainly a question of convenience. 
In this case the choice has fallen upon a classical formula developed by A. R. 
Clarke in his work of 1880 and used extensively in a modified form by J. F. Hayford 
(1909) in geodetic adjustment some thirty years later. 

Referring to Figure 2, the formula defines the meridian component of the 
deflection at the station, produced by a mass of the surface material of the Earth, 
in the form: 


D = 12''-44 x (8/p) x hx (sin 62 —sin 4;) x loge Ro/Ry 


where D is the deflexion caused by a stratum A statute miles thick lying 
within a four-sided compartment of the form shown as a hatched area in 
the diagram. 

5 is the mean density of the stratum 

and p is the mean density of the Earth as a whole. 





Observing 
station 


The value of p in the computations which follow is taken to be 5-515 g/cm? 
and for 5, which in this case must represent the mean density of sea water proper 
to the region, a value of 1-026 g/cm? has been used. 
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Clarke’s formula has been modified slightly in replacing 6 by «, where « is 
measured from south through west so that «a; = 180°—62 and a = 180°—@). Itis 
possible then to define tidal attraction as follows: 


Component of deviation to the south 


1°438 x Cx log, Re/Ri x (sin «:—sin a) x 107% milliseconds 


Component of deviation to the east 
1°438 x C x log, Ro/R) x (cos ag—cos a) x 107% milliseconds 


where C is the amplitude, in centimetres, of the appropriate constituent of the 
tide representative of the compartment. 


By applying these formulae to a network of compartments radiating outwards 
from the station over the adjacent tidal waters up to a range so distant that the 
computed effect becomes of negligible magnitude, it is possible to calculate the 
total component of deflection in the two directions. It is important to note how- 
ever that the phenomenon is periodic and consequently changing in time. The 
time element is introduced, in accordance with normal tidal practice, as a phase 
angle representing the lag of the tidal effect upon the appropriate constituent of 
the equilibrium tide. The phase has a spatial distribution and therefore must 
receive detailed geographical attention. It is necessary to allocate a mean value of 
the phase lag to represent each compartment and to sum the effects given by 
several compartments in a vectorial sense so that the phase as well as the amplitude 
of the total deflection is obtained. 


3. The application of the method 


As an illustration of the method, Figure 3 shows a network of compartments, 
numbered 16 to 33, selected to represent the area of St. George’s Channel, Bristol 
Channel and the South West Approaches, after due consideration has been given 
to the delineation of the coastline and the disposition of the M2 constituent of the 
tide in the area. It can be seen that while the compartments are designed to fit 
the coastline as closely as possible, at the same time consideration is given to the 
variations in the tidal constants. In the north of the area where the variation in Me 
is large, the compartments are correspondingly small while in the south, where 
the tidal changes are small, the compartments in contrast are much larger. 


The calculations are given in Table 1, Section C where the following definitions 
of the column headings apply: 


Column A = the value, allocated to the compartment, of the phase lag upon 
the equilibrium tide of the Mz constituent in the tidal waters. 

Column B = the amplitude of the M2 constituent in the tidal waters proper 
to the compartment and measured in centimetres. 

Column C = the quantity entered under column B but adjusted 
(a) for the inadequacy of the fit of the compartment boundaries 

upon the coastline, or 

(b) for the drying areas within the compartment. 

Ri, Re, «; and a as defined earlier in Figure 2 and in associated comments. 

S = component of deviation to the South in milliseconds. 

E = component of deviation to the East in milliseconds as previously defined. 
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In this manner the Irish Sea, the English Channel, the North Sea and the 
Eastern Atlantic up to a range of approximately 800 miles from Bidston were 
treated as separate regions and both directional components of the respective con- 
tributions of these areas to the total attraction of the mean lunar tide were com- 
puted. It was not considered necessary to proceed beyond this range since the 
contributions from the distant compartments were becoming quite small and were 
further reduced by the introduction of a factor equal to the cosine of the angle of 
depression below the horizon. In any case, in an amphidromic system, such as is 
known to exist to the west in the Atlantic and also to the north near the Faroes, 
contributions from oceanic areas in the middle distance can be regarded as largely 
counterbalanced by the effect of larger tides, in opposite phase, at even greater 
distance. An Atlantic and Arctic limit was selected in this way, partly by experi- 
ment and partly by judgement based upon theoretical argument. 


4- The treatment of local tides 


Clearly the proximity of a large tidal range in the immediate vicinity of the 
station called for special treatment in the computations for the estuaries of the 
Mersey, the Dee and for Liverpool Bay within a range of twenty miles. Unfortu- 
nately it is in this very range that the most serious difficulties are met. On the one 
hand it is necessary to cope with large drying areas in the estuaries where the 
attracting tide appears as a discontinuous feature. On the other hand the phase 
and amplitude of M2 were known precisely only at a small number of locations, in 
fact three in all, and although the three locations provide a good cover in a geo- 
graphical sense, they do not adequately represent the variations in tidal character- 
istics which must be significant in the assessment of attraction. The difficulties 
were eventually surmounted by a detailed investigation of whatever tidal infor- 
mation was available. In the first place differences between high and low water 
times at several localities were converted to M2 phase differences. Secondly, 
although harmonic tidal constituents were not generally available, characteristic 
mean tidal curves were known for certain locations and interpolation provided an 
adequate cover of these characteristic curves. It was found possible to introduce 
sample drying heights into the curves which were then analysed to determine the 
mean lunar tide. This process indicated an effective amplitude and an appropriate 
phase value for Mz under the varying conditions. 

An example of the stages of calculations for a compartment in this troublesome 
region is given below: 


Amplitude of M2 in the compartment 10°oft. 


Drying area as a proportion of the total area of the com- 
partment, obtained by planimeter 0°38 


Mean drying height (above datum) 12-oft. 


On the basis of a 12ft. drying height introduced into the 
characteristic tide curve, itself interpolated to repre- 
sent the compartment, analysis gives the effective 
mean amplitude in the drying area as 72 ft. 


The overall effective amplitude of M2 for the compartment 
is therefore 
(0°38 x 7-2) + (0°62 x 10°0) 8-9 ft. 
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On this basis it was possible to proceed according to the techniques described 
earlier so that satisfactory calculations of the attractive effects for these areas were 
obtained. That this local region should eventually be shown to provide approxi- 
mately half the total effect is justification for the intensive treatment given. 


5. The effect of the curvature of the Earth’s surface 


The depression, below the horizon of the station, of the contributing compart- 
ments was earlier cited as a factor reducing the significance of distant areas. ‘The 
need for a cosine reduction factor in the calculations on this account was a matter 
which required investigation, particularly since the observing station at Bidston 
stands at a height of 178 feet above mean sea level. Fortunately, for tidal waters in 
proximity to the station which have been seen to provide the most important con- 
tributions to the phenomenon under investigation, the reduction effect was found 
to be quite negligible and only on reaching the extreme limits of the area studied 
did these considerations amount to a change of 0-5 per cent in the computed values. 
In view of the relatively small significance of contributions at this range to the total 
effect and in view also of the relative uncertainty in the basic data in these remote 
offshore areas, the reduction factor due to depression below the horizon was 
finally disregarded. 


6. The calculations for the Me constituent 


The calculations have been performed upon a regional basis in Table 1, Sections 
A to F, and the networks of compartments used in the calculations are shown in 
Figure 4 for Section A, Figure 5 for Sections B, C and D, and Figure 6 for Sections 
E and F. For each section the vectorial totals are determined so that the regional 
values of both the magnitude and phase of the deviation can be provided. A final 
vectorial sum of the regional contributions gives the definition of the total vectors 
of the M2 constituent of tidal attraction in Section G. An examination of the results 
shows that the effect of the ocean tide upon the direction of the vertical in any region 
is to introduce a periodic elliptical motion, the delineation of which depends upon 
a number of complex inter-related factors including the azimuth and range of the 
body of water from the observing station, the amplitude of the tide and the phase 
lag of the appropriate constituent in the tide. ‘The axes of the total vector ellipses 
which for Bidston represent the regional contributions to the M2 effect are shown 


in Figure 7 and provide a visual assessment of the attraction in a geographical 
sense. 


7. The inference of the principal solar and larger lunar elliptic semi- 
diurnal constituents, S. and N»2 


Although the disposition of Mz is known with considerable accuracy in the 
waters adjacent to the British Isles, the same cannot be said for any other tidal 
constituent. This is partly due to the fact that whereas it is a relatively simple 
matter to calculate Mz from non-harmonic data, in so far as it represents the 
average tide, it is a much more difficult matter to determine from similar sources 
S2 and Ne, which represent variations upon the average tide. For the latter con- 
stituents one is therefore restricted to a much smaller store of data in the form of 
harmonic constants and even here it is found that the smaller constituents are 
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determined with far less certainty than Me by harmonic analytical processes. For 
these reasons it was considered impossible to proceed for Sz and Ne using the tech- 
niques devised for M2. 

The phase lag and the amplitude of any one tidal constituent varies considerably 
from place to place. It is, however, a well marked characteristic of ocean tides that 
the relationships between the amplitudes and phases of the constituents within a 

















Bidston 





Fic. 6. 


single species show a significant regional consistency. It is because the rather large 
bodies of water in the various channels and gulfs superimpose some restrictions 
upon the nature of the tidal development by virtue of their physical characteristics 
of depth and shape, that the various modifications of the equilibrium tide which 
are found in nature tend to be regional in character. It is on this basis that it was 
found possible to infer the attraction effects of Sz and Nez using mean regional 
values of the amplitude and phase relationships to Mo. 


8. The calculation of the amplitude and phase relationships 


For this purpose, the results of all available harmonic analyses of tidal observa- 
tions in the area were assembled from the files of the Tidal Institute and also from 
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the publications of the International Hydrographic Bureau. The values of the 
amplitude ratios, Mz/S2 and M2/Ne and also of the phase differences (S2—Mg2)° 
and (Mz2—Ne)° were listed for the same six regions which formed the bases of the 
Mz calculations. The task then became an exercise in the systematic reduction of 
data to give the most reliable mean values for each region. 

The procedure was as follows: 


A. The uniformity of the values in each of the regions was examined in order 
to determine whether differences in tidal characteristics were significant 
enough to justify subdivision of the region. 

. If no subdivision was to be made: 
(i) All obviously inconsistent values were first excluded. 

(ii) Weights were applied according to the amplitude of the tide and to the 
length of the tidal record analysed. The weight given was W = M x P, 
where M is the amplitude of Mg to the nearest foot, and where P is the 
number of months duration of the record. 

(iii) Using these weights, average values for the region were determined. 

. Where subdivision was found to be necessary, each subdivision was treated 
as (B) above. The values of the subdivisions were then combined, separ- 
ately for each directional component, weighting the mean according to the 
relative importance of the subdivision in its contribution to the appropriate 
directional component in the M2 calculations. 


An example of the calculations involved is given in Table 2 for the area covered 
by St. George’s Channel, Bristol Channel and the South West Approaches. Here 
it is clear that the Bristol Channel area, taken as subdivision 2, gives characteristic- 
ally different values from those of the rest of the region, particularly in the phase 
relationship of Mz and Sz. Furthermore, the greater coverage of tidal information 
in such a relatively small area would give undue weight to this anomalous feature 
if the region was to be treated as a single unit. In this example, the reduction of the 
data is carried out in full and is self explanatory. Since in assessing the deviation 
of the vertical at Bidston, the Bristol Channel area has less significance in an east— 
west sense than in a north-south sense, then the respective contributions made by 
the two subdivisions to the regional values of the phase and amplitude relationships 
is made to differ accordingly so that a set of values is given for each of the direc- 
tional components. 

In calculations for the North Sea it was necessary to treat the continental 
coastal areas, represented by compartments numbered 1, 3, 5, 8, 9, 10 and 11 of 
Figure 6, as a subdivision giving values distinctive from those of the remainder of 
the North Sea, and when dealing with the English Channel, a more complex 
situation was found which necessitated individual treatment of four separate sub- 
divisions. More uniform characteristics were found in the other three regions so 
that for the Eastern Atlantic, the Irish Sea and for the Mersey, Dee and Liverpool 
Bay, a single set of values can be applied in each case in the calculation of both the 
north-south and the east-west components. 


9. The evaluation of S2; and N2 


The summary of this stage of reduction is given in Table 3 from which it is a 
simple matter to derive appropriate regional results for Sz and Nz based upon the 
Mz values of Table 1. The interpretation of the data in this way is given in ‘Table 4 
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where the regional values have also been summed vectorially to determine the 
total Sz and Ne constituents of deviation experienced at Bidston. Just as the axes 
of the total vector ellipses were produced in Figure 7 for the regional contributions 
to the M2 effect, so Figure 8 reproduces the axes of the total vector ellipses for the 
total Sz and Ne effects and compares them with the corresponding ellipse for Me. 
Finally in Table 5, the fundamental dimensions of all the ellipses determined in 
this study are tabulated. 


Table 4 
Tidal attraction at Bidston 
Mean solar tide So 


To South To East 
Area Amplitude Phase lag Amplitude Phase lag 
millisecs g° millisecs g 

Mersey, Dee, Liverpool Bay 2692 174°5 2°361 I71°5 
Irish Sea 2°177 181-1 2°239 176°0 
St. George’s Channel, Bristol Channel 

and S.W. Approaches 0°621 213°5 0°446 45°6 
English Channel 0166 270°! 0095 357°9 
North Sea 0°365 304°5 0370 171°4 
Eastern Atlantic 0°450 92°2 o-715 351°3 
Total effect 5°186 181-6 3°8907 168-7 


Major lunar elliptic constituent Ne 


To South To East 
Area Amplitude Phase lag Amplitude Phase lag 
millisecs g° millisecs g° 

Mersey, Dee, Liverpool Bay 1°598 108-0 1°401 1050 
Irish Sea 1368 115°7 1°407 110°6 
St. George’s Channel, Bristol Channel 

and S.W. Approaches 0-325 147°2 0238 340°6 
English Channel 0-090 202°4 0-056 288-2 
North Sea O'212 237°0 0°'217 101°9 
Eastern Atlantic 0°250 35°7 0°397 294°8 
Total effect 3°167 115°4 2°440 101°7 


Table 5 


Dimensions of total vector ellipses 


Major semi-axis Minor semi-axis Direction of 
Mean lunar tide, Me millisecs millisecs major axis 
(from North) 


Mersey, Dee, Liverpool Bay 41°-2 West 
Irish Sea 45°°8 West 
St. George’s Channel, Bristol 

Channel and S.W. Approaches 35°°3 East 
English Channel . 2°-3 West 
North Sea 46°-8 East 
Eastern Atlantic 79°°3 East 


Total effect M2 37°°4 West 
Total effect Se . 36°°7 West 
Total effect No 37°°4 West 
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ro. An assessment of the results 


The deviation of the vertical has now been defined as a very small periodic tilt, 
but in order to establish the significance of this phenomenon it is valuable to view 
the results in the light of their implications within the study of earth tides at the 
Bidston station. 


West Approaches 


St Georges Channel 
es Bristol Channel and 


Seuth 


N 


Fic. 7.—Total vector ellipses of the regional contributions to Ms. 


An example can be taken by examining the north-south component of tilt of 
the Me constituent. Observations of tilt recorded by a modified Milne-Shaw 
Seismograph were made at Bidston throughout the International Geophysical 
Year, 1957-8. The harmonic analysis of 19 consecutive sets of 29 day periods of 








‘> East 


Fic. 8.—Total vector ellipses for Mz, S2 and Ne. 


observation gave an observed amplitude of 73 milliseconds for Mz with a phase lag 
of 141°, in the convention that tilts down to the south are positive. This quantity 
is known to contain the body tide of the Earth in addition to those influences intro- 
duced by the oceanic tide. Assuming certain well-established values for the Love 
numbers, it is possible to attribute 5 milliseconds to the body tide at this station 
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with a phase lag approximating to o° and if the body tide is removed from the 
observed quantity there remains a measure of the oceanic influences amounting 
to 77 milliseconds in magnitude and 143° in phase lag. Proceeding with this analy- 
sis, it is known that the oceanic influences must themselves be subdivided into real 
surface tilts due to the loading of the ocean tide in the vicinity of the station on the 
one hand, and apparent tilts caused by periodic deviation of the vertical as a result 
of the attraction of the mass of tidal water on the other. The latter phenomenon is 
here calculated to be of magnitude 17 milliseconds with a phase lag of 138° so that 
the oceanic loading effect must amount to 61 milliseconds with a phase lag of 145°. 
The importance of the effect of tidal attraction in this connexion is thereby estab- 
lished in relation to the associated phenomena and its assessment from a geophysical 
point of view can proceed. 

Already the values given above indicate the lag upon the ocean tide of its loading 
effect operating upon the Earth’s crust. This value can be checked by arguments 
using the east-west component of Me and this, together with the solutions given 
for Sz and Ne are now available for a geophysical interpretation of the continuing 
series of tilt observations at the Bidston station. 


Liverpool Observatory & Tidal Institute, 
Birkenhead. 
1961 March. 
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The Re-magnetization of a Surface Outcrop 
by Lightning Currents 


K. W. T. Graham 


(Received 1961 March 30) 


Summary 


Scattered directions of magnetization are commonly found in 
samples taken from surface outcrops in South Africa. An investigation 
of the directions of magnetization in samples from one outcrop of a 
Pilanesberg dyke combined with laboratory measurements of the 
remanent magnetism produced by d.c. fields and spark discharges 
showed that the pattern of magnetization in the outcrop was due to two 
horizontal electrical currents one of 5 ooo amps, the other of 50 000 amps. 
These currents are almost certainly associated with lightning discharges. 
The effect of a.c. demagnetization in removing the magnetization pro- 
duced by lightning is discussed. 


1. Introduction 


Almost invariably magnetically stable samples from surface outcrops in South 
Africa have shown scattered directions of magnetization while those from deep 
quarries or underground workings have proved to be consistently magnetized. 
(Gough 1956, Graham & Hales 1957, Graham & Hales 1961). Various workers 
have suggested a number of possible explanations for these “‘surface effects” 
including weathering, thermal cycling, magnetostriction, thermo-stress cycling 
and lightning. 

Generally, weathering is detectable and avoidable and the daily thermal 
cycling of an outcrop can be shown to be limited to a very few inches below 
the exposed surface. However, the weathering of the uppermost half-inch of 
rock can produce very high stresses which might be felt several feet down in the 
rock. Similarly the daily thermal expansion of the surface of an outcrop could 
possibly result in a kind of stress cycling. Whether or not these forces could have 
any permanent magnetostrictive effect on the rock is not known. 

Hallimond and Herroun (1933) pointed out that lightning can re-magnetize 
rocks and advised that samples used in laboratory experiments should be taken from 
depths of about soft. If lightning were responsible for the remagnetization of 
rocks, one might expect to find the remanent magnetic directions forming a pattern 
of concentric rings round the point at which the current reached the ground. 
Random sampling over a wide area would be expected to yield samples magnetized 
almost horizontally but with random declinations. The directions of magnetization 
of the numerous samples collected by the author are not concentrated in the 
horizontal plane. It was thought before this investigation was undertaken that it 
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followed that lightning was not of major importance in the remagnetization of 
surface samples. 


2. Work on an outcrop of the Robinson dyke 


The Robinson dyke is a member of the family of Pilanesberg (or Pilansberg) 
dykes. Gough (1956) had shown that this dyke, as well as five others, was randomly 
magnetized at the surface but consistently magnetized at depth. Gough’s figures 
for the basic part of this dyke exposed underground are: 


No. of samples: 42 

Direction of mean N-seeking pole: Azimuth 14°:1, 
Inclination: 71°-8, 

Radius of cone of 99 per cent confidence (x) = 4°:9. 


One of the outcrops of the Robinson dyke sampled by Gough was chosen for 
a study of the “surface effects”. The exposure is about soft square, fairly flat and 
weathered to not more than half an inch below the surface. 

By means of a special drill rig (Grahan & Keiller 1961) a total of 20 oriented 
cores were drilled from this outcrop. Some of the cores were 3 ft long and others 
6ft long. All were drilled vertically. The cores, one inch in diameter, were cut 
into specimens one inch long, the specimens being numbered from the top of the 
core downwards. Generally 10 specimens could be cut from 12 inches of core. 
The direction and intensity of magnetization of each of the 784 specimens were 
measured by means of a spinner-type magnetometer. 


3- Laboratory experiments 


The measurements of the direction of magnetization described in Section 5 
showed qualitatively that the pattern of magnetization could be accounted for in 
terms of horizontal currents near the surface of the rock, and laboratory experi- 
ments on the magnetization of specimens were carried out in order to estimate the 
magnitude of the currents involved. It is convenient to discuss the laboratory 
experiments before describing the results of the field work. 


(a) Magnetization experiments 

Six specimens from the deepest bore cores were selected for experiments using 
a d.c. magnetic field. The specimens had intensities close to 2-8 x 10-$ gauss. 
This value was calculated from Gough’s (1956) figures as the mean intensity of 
underground samples from the basic part of this dyke. 

The direction and intensity of magnetization of each specimen were measured 
before and after each experiment. 

In the d.c. magnetization experiments the specimen was placed in a solenoid, 
the current in which was slowly increased to a predetermined value at which it was 
held constant for a given time and then slowly decreased to zero. ‘The orientation of 
the specimen relative to the axis of the solenoid was in general such that the 
original direction of magnetization was at right angles to the axis of the solenoid. 

The direction and magnitude of the remanent magnetic vector produced by 
the applied field was calculated by vector subtraction from the measured vectors 
of magnetization before and after the application of the field. ‘The known direction 
of the applied field was not used in the calculation. The angle between the direction 
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of applied field and the direction of the calculated remanent vector of magnetiza- 
tion produced by it is a measure of the validity of the method. It was found that 
in all cases this angle was within the limits of the experimental error, never exceed- 
ing 5°. 

This vector of remanent magnetism has been plotted on the vertical axis in 
Figures 1 and 2. It will be seen that the curve starts to depart from the horizontal 
axis at about 8oe and that 120e is sufficient to produce a small but significant 
magnetic vector. The curve shows no signs of approaching saturation at 3800e. 


It would appear as though the best curve through the points has a slight “kink” 
at about 2100e. 
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No magnetic anisotropy was observed in these experiments. 

During the experiments the time of application of the d.c. field was varied 
from 15 seconds to 28 hours. Within experimental error, the remanent vector 
was not affected by the time of application of the d.c. field. The random storage 
of specimens in the Earth’s field for periods of several months showed no significant 
effect. 

In spite of the fact that for this rock the remanent magnetism is independent 
of the duration of the applied d.c. field for periods ranging from 15 seconds 
upwards, it could not be said that the field due to a lightning discharge would 
behave similarly. A lightning flash consists of one or more consecutive strokes, 
usually four in South Africa (Schonland 1956). Each stroke consists essentially 
of a “‘leader’’ (stepped or otherwise) carrying a fairly low current, followed by 
the main flash. It should be noted that although the electrical breakdown for the 
leader proceeds from cloud to ground while that for the main flash is from ground 





88 K. W. T. Graham 


to cloud, nevertheless, the actual current or flow of electrons is in the same direction 
throughout any particular discharge. (Generally the cloud is negatively charged 
and electrons move from cloud to ground.) There can therefore be no alternating 


field as suggested by Rimbert (1958). Instead the field is made up of a number of 
d.c. pulses, all of the same sign. 
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The current in the main stroke varies from 10000amps to more than 160000 
amps. In 50 per cent of the discharges recorded in Rhodesia, the current exceeded 
40000 amps (Anderson & Jenner 1954). The time taken for this current to reach 
its peak varies from 1 to 19s, with an average of 6s, while the time taken for 
the current to fall from its peak value to half its peak value varies from 7 to 115 us 
with an average of 24 us (Chalmers 1957, p. 252). 

With the co-operation of Mr J. Kritzinger of the Department of Electrical 
Engineering of the University of the Witwatersrand it was possible to evaluate 
the effect of very short pulse fields on the magnetic remanence of specimens from 
the Robinson dyke. The rock specimens were placed at known distances from a 
long straight conductor connected via a spark gap to a high voltage impulse 
generator and were subjected to the pulse fields resulting from high current 
discharges along the conductor. Photographs of the trace of a calibrated oscillo- 
graph were taken so that estimates of the current, pulse shape and duration of 
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discharge could be made. (See Kritzinger 1961 for further information on 
circuits.) The direction and intensity of magnetization was measured before and 
after each experiment and the vector produced by the field was estimated in the 
manner described above. The results of these experiments are presented in Table 


1 and the magnetic vector is plotted against the field as large circles in Figures 1 
and 2. 


Table 1 


Specimen Peakcurrent Distance Field Duration* No. Magnetic 
No. amps cm oe ps of shots vector 
gauss x 107% 
5°4 
29 
386 
15°7 
114 
53 
14°5 
13°7 
63 
15°8 
27°5 
33°5 


6/48 4000 16 50 10 
6/48 4000 8 100 10 
6/48 4000 3 267 10 
6/49 4.000 11 73 10 
6/49 4.000 5 160 10 
11/48 3 900 6 10 
11/51 4000 80 10 
11/50 4000 80 10 
11/47 3 200 7° 14 
11/46 3 000 100 1°4 
6/47 3 000 100 1"4 
6/46 3 000 100 1°4 


7 ee ee ee 


* Duration of spark from start, through peak current to half-current. 


The accuracy of this experiment is not as high as that obtainable when using 
a d.c. field generated by a solenoid. The calibration and reading errors involved 
in the estimation of the current should not have exceeded 10 per cent. The error 
in measuring the distance between the centre of the specimen and the conductor 
should not have exceeded 5 per cent but when the specimen is very close to the 
conductor, the curvature of the field and the field gradient across the specimen 
would be appreciable. The specimens were oriented relative to the conductor 
by eye, but the angle between the intended magnetic field and the calculated 
vector reached 9° in one case and otherwise did not exceed 5°. 

It is clear that the effect of d.c. pulses of gus. duration from start, through 
the peak to half current is very similar to that of a d.c. field applied for a much 
longer time. The effect of the shorter pulses of 1-4 us. duration may be a little 
lower but it is not known whether these differences are significant or not. As shown 
in Table 1, the effect of several pulses of the same size is not significantly different 
from that of a single pulse. It is, therefore, very probable that the magnetic 
effect on this rock of one or more main strokes of a lightning discharge will follow 
the curves shown as Figures 1 and 2, and these curves will be used to estimate 
the field and thence the current responsible for the remagnetization of part of this 
outcrop. 


(6) Demagnetization by alternating current 

Alternating current demagnetization experiments were performed on speci- 
mens from this outcrop, primarily in order to determine to what extent, if any, 
it is possible to remove the secondary magnetization without affecting the original 
T.R.M. By this means it might be possible to make an estimate of the direction 
and intensity of the T.R.M. of a rock, in spite of the fact that the original mag- 
netization had been swamped by magnetizations acquired subsequently. 
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The apparatus used is similar in principle to that described by Creer (1959) 
in which the specimen is simultaneously rotated about a horizontal and a vertical 
axis by means of a pair of bevel gears. The alternating magnetic field is slowly 
and smoothly increased to any desired level and then similarly reduced to zero. 
In our apparatus, the lower bevel gear is also rotated but at a speed different 
from that of the main cradle. ‘This ensures a more random motion of the specimen 
which presents every possible direction in it to the a.c. field. The peak field at the 
centre of the coil is set to the desired level by means of a Variac. The coil, mounted 
on a trolley, is then slowly pulled up to and over the spinning specimen by means 
of a winch. After any desired length of time the field at the specimen may be 
reduced very nearly to zero by reversing the winch and removing the trolley to 
about five feet from the specimen. The increase and decrease of the field is, 
of course, not linear with time but it can be made to vary very slowly and smoothly 
over a wide range. It was found empirically that better results were obtained 
using 500 c.p.s. alternating current from a generator than by using 50 cycles from 
the “‘mains’’. A large pair of Helmholtz coils was arranged so as to cancel the 
d.c. field of the Earth at the specimen. It was found experimentally that cancella- 
tion of the Earth’s field or doubling the field made no difference to the a.c. 
washing experiments and the use of the Helmholtz coils has been discontinued. 


N N 


Fic. 3. Fic. 4. 


Lower hemisphere stereographic projections of the directions of remanent 

magnetization of specimens 1/26 and 12/13 after each stage of a.c. 

demagnetization. = = North pole downwards. e = South pole 

downwards. © = Cone of 99 per cent confidence for mean direction 

of magnetization (Gough 1956). Figures denote peak a.c. field in 
oersteds. 


The direction and intensity of magnetization of the specimen is measured 
after each of the successively larger applications of d.c. field. Figure 3, a lower 
hemisphere stereographic plot, shows the way in which the direction of remanent 
magnetization of specimen 1/26 changes with increasing a.c. field. After being 
subjected to a peak field of 1600e its direction of magnetization is very close to 
that determined by Gough (1956) for the basic part of this dyke in Robinson 
Deep Gold Mine. Further increases in the field do not alter the direction appreci- 
ably until the intensity of magnetization is less than about one tenth of Gough’s 
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mean intensity of magnetization for this dyke of 2-8 x 10~% gauss. For larger 
demagnetizing fields, the directions of magnetization start to scatter. 

The initial intensity of magnetization of specimen 1/26 was 17-9 x 107% gauss 
(see also Figure 5) i.e. about six times more strongly magnetized than unaffected 
rock at depth. Similar results are obtainable from specimens with initial magnetiza- 
tions up to about 30 x 10~%, i.e. about ten times the intensity of the original T.R.M. 
In each case the direction of magnetization moves towards and finally pauses at a 
point very close to Gough’s mean direction. Some scatter may occur with fields 
exceeding about 600 0e. 

In specimens whose secondary magnetization is more than about ro times 
stronger than the ‘T.R.M., a.c. washing does not effectively restore the original 
direction of magnetization. Figure 4 shows the progress towards the original 
direction of magnetization of specimen 12/13 whose initial intensity was 2-1 x 
10-1 gauss. No clear “end point’ is reached. After treatment in the a.c. field 
of 3870e, the intensity of magnetization was down to 1°1 x 1074 gauss and it is 
clear that the effective coercive force of such strong magnetization is very close 
to that of the T.R.M. 

Rimbert (1958) drew attention to the difference in the intensity vs a.c. 
demagnetizing field curve for magnetization probably due to lightning and that 
due to a d.c. field. Figure 5 confirms this observation. The demagnetization 
curves for the specimens probably magnetized by natural lightning begin with a 
gentle slope which increases to a point where the curve becomes nearly linear. 


—— Specimen magnetised by Natural Lightning 


Specimen magnetised by D.C. field. 
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—— —Specimen magnetised by Pulse field. 
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The curves for the specimens magnetized by a d.c. field do not show this initial 
gentle slope. Specimens 11/50 and 11/51 were magnetized by the field due to an 
artificial spark discharge (see Section 3a), specimen 11/50 was subjected to one such 
discharge while 11/51 was magnetized by five equal discharges at short intervals. 
The demagnetization curves for these two specimens do not appear to be signifi- 
cantly different from each other and they resemble the curves for the specimens 
magnetized by natural lightning more closely than they do those magnetized by 
a d.c. field. 

The slight plateau formed by each demagnetization curve at intensities ranging 
from 2 to 4 x 1078 gauss probably represents the stage at which the I.R.M. becomes 
comparable with the T.R.M. This supports the assumption that the average 
T.R.M. of this outcrop of the Robinson dyke has much the same intensity of mag- 
netization as the underground exposures studied by Gough (1956). The mean for 
the underground exposures (calculated from Gough’s figures) of 2-8 x 10-3 gauss 
is consistent with these and other demagnetization curves 


7 
1 
Intensity scale (gauss x 10">) 
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Since the addition of magnetic vectors (see Sections 3a and 4b) appears to be 
entirely satisfactory, one would expect that with the removal of one or more 
vectors during demagnetization a similar sort of vector subtraction would be 
possible. In fact this does occur for specimens magnetized by a d.c. field. Figure 
6 shows that at least during the first stages of demagnetization one can picture the 
T.R.M. remaining constant while the I.R.M. is progressively removed without 
changing its direction. The vector resultant decreases in length and changes in 
direction in accord with the parallelogram law and moves along line AB. However, 
Figure 7 shows that in the case of a specimen with I.R.M. due to lightning no such 
simple model is possible. This peculiarity is probably directly related to the initial 
gentle slope of the demagnetization curve for specimens with I.R.M. due to lightning 
but the phenomenon is not yet clearly understood. 


4- The natural remanent magnetization 
(a) The main pattern 

Seventeen of the 20 cores drilled from the outcrop described above came from 
the area about 14ft by 5ft depicted in plan as Figure 8. 

Lower hemisphere stereographic plots of the directions of magnetization of 
the specimens from each of cores 1, 4, 12, 14, 20 and 8 are shown as Figures 10 
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to 15 respectively. On the whole the specimens from these cores as well as those 
from most of the other cores, are magnetized in directions which fall on nearly 
vertical, approximately NE-SW great circles. The “strike” of each of these 
planes of magnetization is shown as the line through each of the “bore holes” 
in Figure 8. As will be seen from the plan, these vertical NE-~SW planes nearly 
coincide with the approximate plane from which the specimens of cores Nos. 8, 
16, 1, 4, 12, 13, 9, 14 and 20 were drilled. 


nan ee Soe 
— —— Initial 1.R.M 


Measured resultant. 
Figures indicate peak A.C. field (oersteds) 
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It is therefore fairly accurate to represent the direction and intensity of mag- 
netization of the specimens as vectors in the plane of a section approximately 
through these cores. Every fourth specimen from each of the above cores is 
represented in such a section which is presented as Figure 9. The magnetic 
vectors for core No. 8, shown as dashed lines, are not direct observations but have 
been calculated as described in Section 4b below. 

Except for core No. 8, the pattern of the direct remanent magnetic observa- 
tions represented in Figure g is qualitiatively consistent with one that would be 
produced by an electrical current flowing along a single straight conductor passing 
perpendicularly into the plane of the section at the point marked C1. However, 
it is equally clear that this simple picture is not quantitatively accurate. For 
example, specimen 16/9 is about 2ft 4in. from Cr and has an intensity of mag- 
netization of 2-2 x 10-2 gauss while specimen 9/17 is very nearly the same distance 
from C1 but has an intensity of 1-9 x 10~! gauss. Many other examples of speci- 
mens to the left (SW) of C1 being much less strongly magnetized than equi- 
distant specimens to the right (NE) of C1 are to be found. 

The simplest explanation for this phenomenon is that there is a gradual 
increase towards the NE in the intensity to which a specimen could be magnetized 
by a given field. This could occur by a progressive change in the quantity or 
nature of the magnetic minerals in the rock. This possibility was investigated 
by subjecting specimens 14/29 and 16/30 to d.c. fields of various strengths. 
Within the limits of experimental error the two specimens behaved identically 
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and in complete harmony with the field vs remanent magnetism curves described 
earlier. 

Two other possible explanations were investigated. One involving the vectorial 
sum of the field due to the current at Cr and another large current in the opposite 
direction but some distance to the right (NE) of C1 was abandoned when a test 
core (No. 20) failed to fit the predicted pattern. The possibility that the current 
at C1, instead of following a thin, wire-like path, spread out in the form of a sheet 
was considered mathematically. It became clear that in any sheet-like model, the 
field at any point below the sheet would be predominantly horizontal even if the 
point is very close to a sheet in which the current distribution is non-uniform. 
Such specimens as 4/5 and 13/5 show that the current path could not have been 
more than one foot wide. 

As pointed out above, the directions of magnetization of the specimens from 
any one core lie on a great circle. The planes of these vertical great circles are 
represented by the lines through the “bore holes” on Figure 8. For any particular 
core the direction of the perpendicular to the current will lie in this plane and its 
sense can be determined by noting the way in which the magnetic vectors rotate 
as one proceeds from the top of the core downwards. This direction for each core 
is indicated in Figure 8 by the arrow heads. It is clear from the non-parallelism 
of these directions that the current did not follow a straight path and a path 
similar to that indicated in the figure seems possible. 

A good approximation to such a curved current path may be made by dividing 
the curve up into a number of straight segments whereupon the problem may be 
tackleu quantitatively in a simple manner. Guided by the intensity of magnetiza- 
tion of some of the specimens nearest to it, a guess may be made at the 
magnitude of the current. The direction and strength of the field due to each 
of the straight current-path segments may be calculated for a point occupied 
by any specimen. With the aid of a stereographic net these vectors in 3 dimen- 
sions may be added. The resultant field can then be converted to I.R.M. by 
means of the curves shown as Figures 1 and 2. The calculated direction of 
the field and intensity of I.R.M. may be compared with the values observed 
for the particular specimen. The current magnitude may be adjusted if 
necessary. 

By trial and error, the current path consisting of the four straight segments 
shown in Figure 8 was found to fit the observations fairly well. ‘The best value of 
the current was found to be 5-0 x 104 amps. 

Table 2 compares the observed direction and intensity of magnetization of 14 
specimens with those calculated on the basis of this model. 

Specimens 5 and 25 of core No. 7 show that the model is not accurate at its 
extremities. Clearly this is due to the assumed semi-infinite straight segments at 
either end of the current path. However, the overall agreement between the 
observed and calculated values of declination, inclination and intensity of mag- 
netization is good enough to indicate to a high degree of probability that the 
pattern of magnetization of this part of the outcrop is due to an electrical current 
of about 50000 amps flowing along an approximately horizontal, narrow but 
curved path. 

The reason for the peculiar directions of magnetization of the uppermost 
three or four specimens from cores 16, 1, 4 and 2 is not known. They may be due 
to very small branch currents which feed the main current. 
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The solid lines in Figure 16 represent the observed directions and intensities 
of magnetization of core No. 8. These directions are also plotted stereographically 
in Figure 15. Note the way in which the direction rotates through nearly 360° 
between the surface and about 3 ft and also the initial decrease in the intensity of 


magnetization to a minimum at about 1 ft gin. (specimen No. 17) followed by an 
increase. 


Table 2 


Observed Calculated 


Specimen D I Intensity D I Field Intensity 
No. deg. deg. gauss x 107% deg. deg. oe. gauss x 10-4 


1/21 243 +12 29 245 +18 103 37 
2/9 244 +55 99 250 +45 137 72 
3/25 256 —5 23°6 261 +10 106 39 
4/13 244 —I 270 253 —I 224 210 
5/13 224 — 36 280 223 —32 261 260 
7/5 215 +1 480 210 +19 662 > 500 
7/25 226 —18 34 224 —7 124 56 
9/13 265 —54 310 275 —5§2 260 260 
12/17 251 —20 146 257 —2I1 195 160 
13/9 256 —48 380 263 —60 377 470 
14/17 282 —50 116 282 —4I1 166 110 
15/17 286 —y 135 283 —7 178 130 
16/13 233 +55 15°2 235 +56 81 21 
20/21 288 —56 41 280 —s1 120 53 
Note: D (declination) is measured clockwise from true North to the North-seeking pole. 
I (inclination) is positive if the North-seeking pole is up. 


If perpendiculars are drawn from the observed directions of magnetization of 
the first 3 specimens, they intersect at the point marked “Current 2”. The 
remanent magnetism of these 3 specimens could be due to a long straight horizontal 
current flowing perpendicularly into the plane of the paper through this point. 
Knowing the intensity of magnetization, the field and thence the magnitude of the 
current may be found. Specimen 1 gives a value of 4-3 x 108 amps while specimens 
2 and 3 give 5-0x10% amps. Specimen 1 may be slightly weathered so that 
5000 amps is probably the better figure. The field and thence the intensity to 
which each of the remaining specimens of core No. 8 would be magnetized may 
now be found. ‘The direction of this vector is given by the perpendicular to the 
line joining the centre of the specimen to the point marked “‘Current 2” and the 
Vector C is shown accordingly in Figure 16 for each specimen. 

On the basis of the demagnetization experiments described above, the thermo- 
remanent magnetization of this core is assumed to be the same as that determined 
by Gough (1956) underground, viz. D = 14°-1 J = 71°-8, mean intensity = 2°8 x 
10-3 gauss and is represented by vector D in the figure. It will be noted that for 
specimens 1, 2 and 3 this vector is negligible compared with the observed vector 
so that the estimate of the current 2 will not be affected by it. 

For any particular specimen, the vectorial sum of C and D does not produce 
the observed vector M, and a further vector A is required to account for the ob- 
servations. This vector, drawn as the dashed lines in Figure 16, may be found by 
completing the parallelogram as shown for specimen No. 8. 

G 





Figure 10, Core 1. 
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Figure 12, Core 12. 
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Figure 14, Core 20. Figure 15, Core 8. 


x = North pole downwards. e = South pole downwards. The 
specimens are numbered from the top to the bottom of each core. 
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The behaviour of vector A is rather more systematic than that of the observed 
vector M. Its direction remains more or less constant and its intensity falls off 
smoothly with depth. Figure 9, in which this calculated vector A is shown by 
dashed lines, shows clearly that it fits in with the rest of the “Main Pattern’”’. 
It is undoubtedly due to current 1 described above. 


5. Discussion 


Both the electrical currents responsible for the remagnetization of part of this 
outcrop followed predominantly horizontal paths close to the surface of the rock. 
The magnitude of the currents, their behaviour and the characteristic a.c. 
demagnetization curves leave very little doubt that they are associated with 
lightning discharges. 

From Figure 8 it will be seen that current 2 could possibly curve towards the 
North to join up with the main current. Another small current may have joined 
the main current near core No. 15. This could possibly explain the magnetism 
of core No. 6. It seems probable that cores Nos. 10 and 11 were magnetized by 
another large current to the West of them. This current may or may not have been 
associated with the main current described above. 

In the area of the outcrop studied, the effect of a vertical lightning current is 
not seen at all. It seems probable that the horizontal currents could remagnetize 
a greater area of rock than would the vertical current with which they are associated. 
It is clearly fallacious to argue that if the pattern of magnetization of a surface 
outcrop is not one of concentric horizontal vectors, then it could not have been 
caused by lightning. 

For this rock the minimum field required to produce a significant change in the 
remanent magnetism is about 120e. It can be shown arithmetically that to avoid 
entirely such fields due to lightning one should not use samples from a depth of 
less than soft. Securing accurately oriented samples from a depth of soft is a 
formidable task. Deep quarries, road cuttings and mines should be used wherever 
possible but in South Africa such workings do not often exist in the rock formations 
of most interest to palacomagnetism. 

The next best technique at present available is that of drilling to 3 or 6 feet. 
Although the directions of magnetization at these depths can be very significantly 
different from that of the original N.R.M., a.c. “‘washing”’ would, in most cases, 
remove the unwanted effect. 

Clearly, if either vertical or horizontal lightning currents have affected the 
remanent magnetism of an outcrop, systematic sampling over an area of say 50 ft by 
50 ft cannot yield a consistent, but wrong impression of the original magnetic 
direction. The results will merely be negative in that they will be scattered. 

Having shown that the overall random magnetization of one outcrop is due to 
lightning currents, the question arises whether this kind of remagnetization is 
common or not. Gough (1956) sampled ten outcrops of five Pilanesberg dykes, 
including the outcrop sampled by the present author. All ten were found to have 
scattered directions of magnetization. There is no reason to believe that the cause 
of the scatter was anything other than lightning currents. The marked decrease 
in scatter of the underground Karroo dolerites when compared with the surface 
samples (Graham & Hales 1957) has already been mentioned. Similarly samples 
of the lower “‘Shales” of the Table Mountain Series taken from a deep road cutting 
were found to be consistently magnetized whereas those from the mountain side 
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were inconsistently magnetized (Graham & Hales 1961). In both these cases 
it is thought that lightning currents may have been the cause of the scattering. 
Workers in other countries have frequently encountered surface outcrops with 
scattered magnetization. These negative results are generally mentioned only in 
passing, or are not reported at all so that it is impossible to form an opinion as to 
the frequency or cause of such phenomena. However, the “‘minor deviations” in 
the magnetic directions of certain outcrops of a family of Tertiary tholeitic dykes 
of North England reported by Bruckshaw and Robertson (1949) may be attribut- 
able to lightning currents. The patterns of magnetization of some outcrops studied 
by Bruckshaw and Robertson are very similar in character to that described above. 
The possibility of the remagnetization of many surface outcrops by lightning 
currents makes it necessary to use the “conglomerate stability test” with caution. 
J. W. Graham (1949) argued that if each boulder or rock slab in a conglomerate 
is consistently magnetized but the remanent magnetic directions for the various 
boulders is different to the point of randomness, then the rock from which the 
boulders are derived is probably magnetically stable. However, it is obvious that, 
particularly where horizontal lightning currents are involved, the direction of 
magnetization of a surface outcrop may vary rapidly over a small area. Thus, this 
test is only valid if either (a) the samples came from a considerable depth below 
the surface or (b) it can be shown that the separation of the specimens from each 
boulder is of the same order as the distance between the various boulders sampled. 


6. Conclusions 


The magnetic fields generated by large electrical currents associated with 
lightning discharges have re-magnetized the portion of the outcrop studied and 
may be the cause of many of the “scattered” directions of magnetization encountered 
in the palaeomagnetic study of surface samples. The only completely satisfactory 
way of avoiding this hazard is to restrict oneself to samples collected from depths 
of greater than about soft. However, systematic sampling over even a fairly small 
area cannot yield consistent but “wrong”’ directions of magnetization resulting from 
lightning currents. If scattering is found, a.c. washing may prove useful, parti- 
cularly if the specimens have come from moderate depths such as 3 to 6 feet. 

J. W. Graham’s (1949) “‘conglomerate test’’ for magnetic stability should be 
used only with extreme caution. 
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Note added in proof.—A recent paper by A. Cox (U.S. Geol. Survey Bull. No. 
1083E, 1961) further demonstrates that the re-magnetization of surface outcrops 
by lightning currents is very common. 


Bernard Price Institute of Geophysical Research, 
University of the Witwatersrand, 
Johannesburg, 
South Africa. 
1961 March. 
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On Solar Wind Generation of Polar Geomagnetic 
Disturbance 


K. D. Cole 


(Received 1961 March 24) 


Summary 
Solar wind generation of polar geomagnetic storms is examined on a 
model suggested by Piddington. His theory is criticized and in- 
ferences are drawn which may be tested against observation. A com- 
parison of solar wind and atmospheric dynamo generation of magnetic 
disturbance (and auroras) is made. It is suggested that a combination 
of these two processes may be necessary to explain the total phenomena. 








Fic. 1. (Taken from Fukushima & Oguti 1952.) Current-system of the 
Sp-field for the equinoctial season. (Electric current of 1-5 x 10° amp 
flows between adjacent stream lines.) 


1. Introduction 


Piddington (1960) uses the hypothesis of a “solar wind” blowing outwards 
from the Sun in his explanation of the longitude dependent geomagnetic distur- 
bance Ds (cf. Figure 1). (A related suggestion comes from Johnson (1960)). It is the 
purpose of this paper to re-analyse his model with the view to criticizing his analysis 
and to drawing inferences from the theory which may be tested by observation. It 
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will be assumed that the reader is familiar with Piddington’s (1960) version of 
analysis of the model. 


2. General description of physical processes 


It is assumed that the “solar wind’’ does not invade a cavity (whose surface is 
S) about the Earth. The wind blows with a velocity w, through the geomagnetic 
field H. On Piddington’s (1960) theory one would expect S to be elongated away 
from the Sun and dome-shaped on the day side. Figure 2 illustrates a cross-section 
of S by a plane normal to the Sun—Earth line through the centre of the Earth. 
Such a cross-section may be about 5 Earth radii in radius. 


Fic. 2. Schematic representation of cross-section at right angles to 
Sun-Earth line of the cavity S formed by the solar{wind in the Earth 
environment. 


Let us consider the instant t = o at which the solar wind has just commenced 
blowing (see Figure 2). The motional induction field in the solar wind E;{ = w a H) 
causes a charge distribution over S. The polarization field of this charge has a 
value E, in the wind and a value E,’ inside S. Ep opposes E; but does not neces- 
sarily balance it. E,’ is unopposed and, having a component along the geomagnetic 
field, causes current to flow along it. This explains Piddington’s (1960) equation 8 
for which he gives no explanation; this gives reason for the j(= 0, 0, j.(r)) of his 
section 6 which he assumed exists (as a twist of the magnetic field). It will be seen 
later in Section III(a) how this electric field Ey’ is essential to the phase propagation 
of such a twist. 

Between S and the ionosphere (I) the current flows principally along H for the 
ratio of the conductivity along and across H is extremely large (~ 104 or 10°). The 
flow of current along H is terminated in the ionosphere where the consequent 
build-up of space charge causes current to flow across the magnetic field. (In the 
whole of the region inside S the best conducting path across the magnetic field is in 
the dynamo region of the upper atmosphere ~ 100 km height above the Earth’s 
surface.) The magnetic field lines between S and the ionosphere are then vir- 
tually electrostatic equipotential lines. 

In the initial stages of establishing a steady electrostatic field system in the 
ionosphere the phase speed of current flow along H inside S is the Alfvén speed as 
demonstrated by Piddington (cf. Section 3(a)). In the steady state it is the electro- 
static field system established within the ionosphere which drives current there 
responsible for geomagnetic variations as observed at the ground (cf. Section 


3(b)). 
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3- Theory 
(a) Current flow between S and I 


The current density (j) in a fully ionized gas in which the gradient of pressure 
of electrons can be neglected is given by (Cowling 1957) 


ooE’ = j+ wrH“ ja H, (1) 


E,’+vAH 

electrostatic field 

momentum per unit mass produced inside S. 
gyro frequency of electrons 

time between collisions of electrons with ions 
ne=r/m, 

charge, mass of an electron 

electron density 

magnetic field. 


j E,’ E,’ — 
= _+ \+ ’ 
J a0 OE, +02 


o| = ao/(1 + w*7?), 


o2 = GWT, 


and ||, | represent components parallel, perpendicular to H. 

In the Earth’s exosphere between S and J, w may vary between 10? (at J) and 
104 (at S) whilst + may vary between 107% (just above J) and 1 (at S). Thus be- 
tween S and J, wr > 10%, thus o9 > o2 > oa). Hence equation (2) is well approxi- 
mated by 


: = eam? 
j = %E,'+o2——-, (3) 


H 
and the equation of motion 


m H 
a PR im 


where p is the density, becomes, using equation (3), 


dv , 
PR = oHE,’. (5) 


Now consider the following model. In a system of cylindrical coordinates 
(r, 0, z) a magnetic field H = (0, 0, H;) exists prior to the application of an electric 
field E, = (0, 0, E,). The field is assumed to contain a uniform distribution of 
fully ionized gas. As a result of E, currents will flow causing additional magnetic 


field and motion. At this subsequent time let v = (v;, vy, vz) and H = (Ay, Ay, 
H,). 
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Then equation (5) yields 


dv 

= = —ohE,+ o2H{v,H.—vzhy), 
dv, 

oy = ogh,Ez+ o2H {vzh,— vyH2), 
dv 

— = = o2H(vrhy— Uf). 

To the writer’s knowledge the equations (6), so fundamental to present-day 
research in controlled nuclear fusion, have never been solved completely. It may be 
possible to get an inkling of the morphology of the processes they describe in certain 
particular cases with the assistance of physical intuition rather than mathematical 
deduction. 


In the idealized twist propagation examined by Piddington a tacit assumption 
has been made that v,; = hy = 0, viz., there is no pinch effect. In this case the 
equations (6) yield 


Vplhy = v2/Hz, (7) 
and 


dv, oh 8 
P dt oh Ez. (8) 
In equation (7) vz appears to be arbitrary. However, it represents the relative 
velocity of matter and magnetic field and therefore is the phase speed (i.e. Alfvén 
speed) of the (vy, 4,) disturbance along Hz. It is clear that the j (= 0, 0, j.(r)) of 
Piddington’s section 6 is due to the component of electric field along H, (cf. 
equation 3) so that j,(r) = ooF,. Piddington offers no explanation of this current. 
Moreover it is seen from equation (8) that the current (providing a j a H force) 
which changes v, is the Hall (c2) current generated by E, acting across h,; this 
explains Piddington’s J,. Piddington’s explanation of his J, is confusing, first 
because he invokes isotropic conductivity and secondly because his equation (15) 
is inconsistent with his equation (13) and the statement to which it is connected. 
Piddington’s twist corresponds to a special case of equations (7) and (8) above in 
which v, takes a stepwise change across the wave front (a phase surface). 
Piddington suggests that independent “twist cylinders” develop in the dawn 
and dusk longitudes and “grow” (at the Alfvén speed) towards the Earth along the 
geomagnetic field lines. The above analysis shows that a more reasonable sugges- 
tion is as follows: surfaces (Spn) of particular phase (in disturbance) slough from 
S and shrink towards the Earth at the Alfvén speed along the geomagnetic field 
lines. This means that the Hall current suggested by equation (8) and likewise 
Piddington’s J, current find continuity on these surfaces (Spn). This dispenses 
with the necessity of having a return current back along the surface of the twist 
cylinders suggested by Piddington. It also nullifies his explanation of the pre- 
liminary reverse impulse Ds(PR) which he attributes to a polarization field (his 
Ey») due to partial inhibition of J,. For, since J, is continuous over an Spn, this 
Ey no longer exists. A more likely explanation of Ds(PR) (which is inevitable on 
Piddington’s theory) is that the sign of E,’ initially established in the ionosphere 
may be sometimes opposite to that of the finally established field. This would 
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depend on the distribution of ionization inside S at t = o. This would explain 
also its generally small size. 


(b) Electrostatic field in the ionosphere 


From t = o onwards current will flow in the ionosphere due to E,’. The 
steady state will be reached only after a time equal to the transit time of a hydro- 
magnetic wave from S to J. However, since the instant ¢ = 0 is indefinite to the 
extent of the time it takes the “‘solar wind” to traverse the system, there can be 
little point here in distinguishing these events except in so far as it bears on the 
theory of Ds(PR). It is the opinion of the writer that the propagation of phase 
discussed in Section III(a) above is substantially irrelevant to enlightenment on 
the direct cause of the main features of that Ds current system which results from 
the model considered here. This will now be discussed. 

In the steady state the geomagnetic field lines from S to J will be virtually 
electrostatic equipotential lines. Current will flow orthogonal to the magnetic 
field lines but in greatest concentration in the lower ionosphere (cf. Baker & 
Martyn 1953). Thus current generated at S will flow down the geomagnetic 
field lines in the dawn hemisphere (under the influence of a very small electro- 
static field), across the ionosphere principally in polar regions then up the geo- 
magnetic field lines in the dusk hemisphere. There may also be some ionosphere 
current at lower latitudes with continuity either in the same or opposite hemispheres. 

In order to determine approximately the morphology of the electrostatic field 
established in the ionosphere we first assume that the space inside S is an insulator. 
This enables the charge distribution on S to be readily determined in a steady state. 
Allowance is then made for the fact that the space inside S is a good conductor. 
Thus the charge distribution that would have formed on S diffuses along the geo- 
magnetic field lines under the influence of an external electrostatic field. The 
sources and sinks of electrostatic field in the ionosphere may then be found 
approximately. 

The current flowing in the solar wind is 


, HaE, 
j= ooE, + mE, + 02 H A (9) 


where E = E,+waH and Ey is the polarization field in the wind. In the steady 
state (/@t = 0) no current would flow at the surface S in the direction of the unit 
normal (n) to S and Ey would be entirely normal to S, viz.j - m = o. Thus on the 
assumption that S is an insulator we find using equation (g) that at S 


-mHawaH -n/H+owaH-n 
~ Go(H » n)2/H2+ 0;(1 —(H - n)?/H2) (10) 
Since oo > o; here, we can discuss equations (10) for two kinds of regions: 
(i)H-n wo, 





Ep 


—-oHawaH-n/H+o;waH-n 





Ey = 


o} 
(ii) H-n #0, 
—oHawaH- n/H+o;waH-n 


"aS oof - n)?/H? 
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It follows from equations (11) and (12) that since o9 > o; the most dense accu- 
mulation of space-charge on S occurs around points such that H- n = 0, i.e. 
where a line of force is tangential to S. With the shape for S suggested above we 
should find two lines along which H+ n = o: One line in the dawn hemisphere 
would possess positive charge, and the other line in the dusk hemisphere would 
carry negative charge. 

Now the space inside S is not an insulator but consists of a system of highly 
conducting electrodes (the geomagnetic field lines). As a first approximation one 
would expect the electrostatic field created in the ionosphere (in the steady state) 
to be that obtained by projecting the charge distribution on S along the magnetic 
field lines to the ionosphere. Hence there will be a source of positive charge at a 
line in the dawn ionosphere such that the geomagnetic field lines through it are 
tangential to S. Likewise there will be a sink of charge in the dusk ionosphere. 
Thus instead of just two relatively small areas in each polar ionosphere for the 
precipitation of opposite space charge (Piddington 1960) there would be two areas 
elongated in longitude and thin in latitude extent. 


(c) Current flow in a uniform, horizontally stratified ionosphere. 


The electrostatic field patterns in the polar ionosphere may be obtained approxi- 
mately by regarding this region as plane with an homogeneous orthogonal (vertical) 
magnetic field. The electrostatic field lines will be coplanar with the ionosphere. 
In the quiet ionosphere the ratio (y) of the height integrated Hall and Pedersen 
conductivities, viz. {a2 dh: {o; dh is about 2 (c.f. Piddington 1960). Thus by equa- 
tion (2) the current direction makes an angle of tan~!) with the electric field 
direction. 

Piddington suggests that the current from a radial electric field would flow in 
closed loops (c.f. the paragraph containing his equation 2). Rather the current flow 
due to a radial electric field would be along centred equiangular spirals r = 
a exp(@/n). Now from the known structure of the ionosphere Hall current pre- 
dominates over Pedersen current around 100 km height but Pedersen current pre- 
dominates above about 130km (cf. Baker & Martyn 1953). ‘To express the 
magnetic effects at the ground in terms of an equivalent current system (at a 
specified height) we must weight the Pedersen current by a factor on account of its 
greater height. For the quiet ionosphere assumed by Piddington an underestimate 
of the factor would be }. Hence the equation of the spirals of the consequent 
equivalent current system would be r = a exp(6/4). Such spirals increase their 
radius vectors by a factor of exp(27/4), i.e. x4 in one turn. Such curves are 
vastly different from the closed loops Piddington attempts to explain. Additional 
assumptions about the structure of the ionosphere would have to be made to make 
the spirals approximate loops, for instance, that the ratio of total ionization 
below 150 km to that in the F region is very large. An idea of a lower limit to the 
ratio 7 in the equivalent current system can be got as follows. There are approxi- 
mately 30 stations spread over an area of 3 000 km by 3 000 km to determine the 
structure of the Ds current loops. This gives an average of about 500 km between 
stations. One would therefore expect a resolution of 500/3 000 = 0°17 in deter- 
mining the dimensions of loops. So that to approximate loops spirals would be 
allowable with an opening factor of 1-17. Thus a lower limit of » = 27/In 1-17 ~ 40 


is established. This is an order of magnitude greater than in the quiet ionosphere 
assumed by Piddington. 
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It has been observed above (Section 3(b)) that the source and sink of charge 
in the ionosphere will be elongated in longitude and thin in latitude extent. How- 
ever, their precise elongation depends not only on the shape of S but also on the 
distribution of v. There appears little hope at present of finding either of these 
two things. Hence to a first approximation the source and sink of charge are con- 
sidered as points on the dawn and dusk longitudes respectively. 

Considering the polar regions to be flat with a vertical magnetic field we may 
approximate the electrostatic field by that of a dipole » ~} for the equivalent 
current system in a quiet ionosphere (assumed by Piddington). Figure 3 shows the 
two systems of interlocking spirals of current that result from a dipole electric 
field. These are vastly different from the closed loops that Piddington wished to 


explain. 
og 


Fic. 3. Flow lines of current system that would be produced by a 
dipole electric field caused by a solar wind. (Magnetic field into plane of 
paper.) 








a 








It may be that in fact the Ds current system is an interlocking pair of systems 
of spirals (cf. Figure 3) and that the present network of observing stations is not 
good enough to resolve them. It may be of significance that spiral forms of the 
aurora do occur. The above theory suggests that such spiral auroras should open 
clockwise in the southern hemisphere (as seen from the ground) and anti-clockwise 
in the northern. This should be checked against observations. 


4- Comparison of solar wind generation theory and dynamo theory of 
aurora and magnetic disturbance 


On the basis of observations then available Cole (1959, 1960) suggested that 
auroras constitute luminous segments of lines of the Ds current system (cf. Figure 1). 
Observations reported since (Davis 1960; Denholm & Bond 1960) support this 
view. In the following discussion no distinction in plan position will be made 
between aurora and disturbance current; they will be treated as aspects of the one 
phenomenon. A simple model of the pattern of auroral orientation and movement 
that is not inconsistent with the reported observations is as follows: closed loops of 
the Ds system contract about a morning “cell” and expand about an evening “‘cell”’. 
(Since the observations of Davis & Kimball 1960, were taken from near the lati- 
tude of the centres of the cells they may be explained by fluctuations in latitude 
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of the cells’ centres.) This approach gives primacy in physical importance to the 
current system. Luminosity (in the form of auroras) may be brought about in at 
least two ways. Firstly, a small variation in the Ds current (brought about by a 
change in conductivity or electric fields in the ionosphere) may cause a hydro- 
magnetic wave to propagate upwards into the region of geomagnetically trapped 
radiation and scatter protons downward (cf. Welch & Whittaker 1960) thus 
producing luminosity and further magnetic variation and so on till there is deple- 
tion of “available” trapped radiation. Secondly, the Ds system may act simply as a 
large magnetic anomaly comparable in size to the South Atlantic anomaly, causing 
the lowering of mirror points of trapped charged particles into the atmosphere 
with similar consequences. On this physical basis the auroral isochasms and the 
diurnal variation of auroras may be explained (cf. Cole & Jacka 1961). 

It is of interest to compare the above electrostatic field theory and Piddington’s 
(1960) theory with a dynamo theory of the aurora and magnetic disturbance 
developed elsewhere (Cole 1959, 1960). 

Piddington (private communication) states that it is the induction field (his E;) 
which drives Hall current in the ionosphere and not space-charge accumulations as 
he stated previously. This will be explained more fully by him in a sequel paper. 
However we can discuss the effects of an electric field in the ionosphere without 
being concerned with its origins (whether it is the E; of Piddington or the E,’ 
deduced above). 

We have already seen that an electric field of roughly radial configuration 
cannot explain closed loops of the Ds system in a uniform ionosphere unless there 
is a great enhancement of ionization in the E region. (There is no conclusive evi- 
dence that this condition of the polar ionosphere is as general as magnetic dis- 
turbance.) On the contrary the polar ionosphere is far from uniform as postu- 
lated in Section 3(c) above and by Piddington. The evidence is strong that some 
auroras contain abnormally high current density compared to their surroundings 
(cf. Heppner 1954). Auroras frequently appear as striations on ionograms. Radar 
observations of some auroras is common. Rocket investigations of the polar 
ionosphere show large gradients of electron density in the vicinity of auroras. 
These things support the assertion that the Ds aurora system has considerable 
filamentary structure in the form of numerous concentric and successively larger 
loops of high line density of ionization (cf. the lines of Figure 1) superimposed on 
a background of diffuse ionization. The current created in such structures by 
electric fields is vastly different from that in a uniform ionosphere. 

Let us consider then an auroral form (say a homogeneous arc) as an infinite 
plane slab of highly ionized air immersed in the ionosphere parallel to the (vertical) 
geomagnetic field. An electric field is applied perpendicular to H (Figure 4). It 
has been shown elsewhere (Cole 1960) that the component of velocity (V) of the 
slab in the direction of the unit normal (m) to the slab is given by 


V-n z= E®%sin a/(1 + ver4/wewi), (13) 


where w,; are the gyro-frequencies of electrons, ions and v,,; their collision 
frequencies with neutral particles. The current (jp’) flowing along the model aurora 
is given by equation 8(b) of Cole (1960). In the steady state and when the electron 
density in the aurora is very much greater than that outside it this equation yields 


jp = —o3E° sing (14) 


in the direction of n a H, where o3 = 6) + 622/03. 
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Figure 5 shows plots of the factors 1/(1 + vevj/wew;) and og as a function of 
we/ ve and height assuming with Baker & Martyn (1953) that (we/ve)(w4/v;) = 1 650. 
Equation (14) demonstrates that an aurora moves normal to itself (i.e. laterally) 
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Fic. 4. Schematic representation of an aurora as a plane slab of highly 

ionized gas (between the parallel planes) parallel to the geomagnetic 

field H in the ionosphere. An electric field E°® is applied in direction 
making angle « with unit normal n to the auvora. 
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Fic. 5. Approximate plot of (a) o3/osmax in an auroral model with 

constant electron density at all heights; (b) the factor 1/(1 + vevs/wews) as 

functions of we/ve assuming with Baker and Martyn (1953) that 
(we/ve)(ai/ va) = I 650. 


provided the electric field has a component (horizontally) along the aurora, whilst 
equation (15) shows that current flows horizontally along the aurora providing 
there is a component of electric field in that direction. 
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Let us represent one of the Ds loops as the circumference of a circle with high 
line density of ionization. If a radial electric field is applied (as in the case of sup- 
posed solar wind generation) then in the steady state there will be neither current 
along nor movement of the loop (by equations 13, 14 above). This occurs because 
the Pedersen current is inhibited by polarization and when this occurs all current 
flow ceases. 

According to dynamo theory it may be suggested that radial flow of wind 
would occur thus causing azimuthal electric field. Thus in the steady state there 
would be radial movement of the loop (equation 13) and current flow around the 
loop. This occurs because the loop causes polarization of the Hall current which in 
turn builds up the circumferential current. Thus dynamo theory seems better 
able to explain the observations. 

Let us modify the strictly circular loops just discussed into the (filamentary) 
loops of a Ds (aurora) system (see Figure 1). Then from solar generation one expects 
a source of electric field at the centre of the morning cells and a sink at the centre of 
the evening cells. Thus the pattern of auroral movement to be expected in a frame 
of reference fixed with respect to the Sun is schematically as shown in Figure 6(a). 


o” 





42" 12 
(a) (b) 
Fic. 6.—Schematic representation of pattern of lateral auroral move- 


ment expected (a) from solar wind generation theory; (b) from the 
author’s atmospheric dynamo theory. 


This pattern is certainly contrary to the observations (cf. Bond 1960) of the 
equatorwards lateral movement of auroras prior to midnight outside the auroral 
zone. On the other hand the dynamo theory would postulate a source of wind at 
the centre of the evening Ds cell and a sink at the centre of the morning cells. The 
expected pattern of lateral auroral movement would then be as schematically 
presented in Figure 6(b) (ina frame of reference fixed with respect to the Sun). 
Observations available so far do not contradict this scheme. 

A difficulty with an ionospheric electric field of the type involved in Section 
3(b) above is that of explaining the uniformity of auroral movement at all heights. 
For by equations (13) and (14), since E® would be practically constant at all heights, 
both V- mand jp’ would vary considerably with height (cf. Figure 4). This difficulty 
is avoided in dynamo theory (Cole 1960) for there the operative electric field and 
movement is determined by a limited height range. 
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Dynamo theory is not without its difficulties. Piddington (1960) has enumerated 
some of these. Basic difficulties are explanation of the mechanism of establishment 
of the necessary upper atmosphere wind system and its rotation with the Sun, 
and of the necessary patterns of ionization. 

It may well be that the ultimate explanation of the morphology of Ds and 
auroras could involve both dynamo theory and electrostatic field theory. Certainly 
both appear capable of explaining magnetic effects of similar orders of magnitude. 
Winds of adequate magnitude in the atmospheric dynamo region are certainly 
observed (Elford 1959) and some electrodynamic consequences of the type analysed 
by Cole (1960) must be expected. On the other hand there are features of the 
electrostatic field theory which are highly attractive, viz. its explanation of de- 
pendence on the position of the Sun and its ability to explain discontinuous changes 
in magnetic disturbance. But it has its roots in an, as yet, unobserved solar wind 
the existence of which is not proved beyond doubt. 


5- Discussion 

Elsewhere the writer (Cole 1961) has cast doubt on the existence of a “‘solar 
wind”’ interpreted in the theory of Parker (1957) as an expansion of the solar corona. 
Rather the writer has suggested that the inner corona is a hydromagnetic wave 
generator and that the expulsion of such waves would tend to inhibit its expansion. 
However, the above theory of “solar wind’’ generation of geomagnetic disturbance 
is not at variance with the concept of hydromagnetic wave emission by the Sun. 
For, assuming an Alfvén speed of about 10% cm/s in interplanetary space, waves 
of periods 2 100s would have wavelengths = 10 Earth diameter. All periods of 
waves up to a day or more may be represented in the emission spectrum of the 
inner corona. Such waves would cause wind of variable direction about the Earth 
in interplanetary space. Such variable winds may account for the complexity of 
auroral movement and magnetic disturbance. Perhaps the strongest objection that 
can be made to solar wind generation of magnetic disturbance as developed by 
Piddington (1960) and in the analysis above is that the “average” Ds system which 
it attempts to explain is the average of a large number of elemental disturbance 
systems (Nagata & Fukushima 1952) which differ significantly in structure from the 
average. A physical theory must explain both the element and its average. This 
seems within the bounds of atmospheric dynamo theory. 

The need for precise determination of the instantaneous morphology of the Ds 
current system and of the simultaneous global pattern of auroras is vital to the 
further development of theory. 
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Small Corrections in the Theory of Surface Waves 


Harold Jeffreys 
(Received 1961 April 13) 


Summary 


Rayleigh’s principle is used to derive expressions for the group 
velocities of Love and Rayleigh waves without numerical differentiation. 
It is also used to find expressions, accurate to the first order, for the 


effects of small changes of the elastic properties on the wave and group 
velocities. 


1. Introduction 


The standard method of determining group velocities of surface waves is to 
assume values of the wave-velocity c and of «x, where 27/« is the wavelength; 
solve the differential equations; and for each « or c interpolate a value of the other 


so that the boundary conditions are satisfied. Then the group-velocity C is 
determined from 


C d d / dx 
oe ae ae 


One disadvantage of this procedure is that for any new suggested structure 
the whole process must be carried out ab initio, becoming more complex with each 
additional complication suggested. Another disadvantage is that the differentiation 
to give C has to be carried out numerically. But the whole range of c is small, and 
the intervals used cannot be less than o-1 of the range of c to give an adequate 
determination of the behaviour of C; in practice two figures are lost in the differ- 
entiation. 

E. Meissner (1926; see also Jeffreys 1959, p. 110) gave a method for Love 
waves, based on Rayleigh’s principle, that avoids numerical differentiation. This 
method could be adapted to Rayleigh waves, but this does not appear to have been 
done. Rayleigh’s principle can also be used to estimate the effects of small changes 
of the distributions of density and elastic moduli. Since it is fairly easy to work out 
the solution for a model that can be regarded as a useful first approximation, it 


appears that such a method should permit rapid testing of any reasonable suggested 
structure. 


2. Love waves 
If the displacement is 


Uj, U2, Ug = O, (2) cos(yt— Kx), Oo 
115 ; 
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we have, apart from constant factors, 


aT = ( py2e? de, 2V = | p(x?2v2 + v'2) dz, 


and by equating 7 and V, 


Ioy? = Nx? +s (3) 
where 


Ih = | pv? dz, I, = | words, r seer” z. (4) 


Then Rayleigh’s principle says that if we work out (4) using a slightly incorrect 
form of v, then for given « (3) will give y? with an error of the second order. 

If we change xo to xo +4« but use for x = xo+8« the form of v calculated 
for xo, this has a first order error, but 


lo(yo+ 5y)* = Lh(Ko + 8x)? + Ip + O(8x)? (5) 
whence to the first order in 6x 
Igyody = Tyxodk. (6) 
But since u/« = c, dy/dx = C, this is 
cC = Ih/h. (7) 


This is Meissner’s result, which is exact. It replaces numerical differentiation by 
the integration of two known functions. 


But in the same way, if we also vary p and yu, and retain the same values of v, 
we shall have 


(Lo + 51o)(yo + Sy)? = (i + 51)(ko + 5x)? + Le + Sle + O((8p)?, (du)?, (5x)*) (8) 


where S/o, 54), 5/g are independent of x. Hence if 5« = 0, 


(Lo + 51o)(yo + by)? =2 (h +8]))«o? + Io + 512, (9) 


determining c with a second-order error for a fixed «; and from (7) 


(co+ 8ce)(Co + 8C) = (1 +6h)/(lo+ 6h), (10) 


giving C also with an error of the second order. 
A method of mine (Jeffreys 1928) used Rayleigh’s principle to allow for increase 


of the velocity of S with depth in the lower layer, by means of (5); but C was then 
determined as usual by numerical differentiation. 


3- Rayleigh waves 
The displacement is 


Uj, U2, Ug = u(z) sin(yt—«x), 0, w(z) cos(yt— Kx). 
We find 


aT = | py"(u? + w?) dz 


2V = [ {A(eo’ — xu)? + 2y(m2u? + w'2) + po(u’ + ew)?} dz 





Small corrections in the theory of surface waves 


ylo = «2I+2xle+Ts 


Ip = | plu®+ w®) ds 


h= [{Q+ 2p)u2 + p02} de 


I, = | (—Aw'u+ pu’w) dz 





Ig = | {(A+2)w’? + pu'?} do. 
Varying « as for Love waves we have 


IgcC = 1, + Ie/x. 
If we vary p, A and yu, without varying « or the trial solution, we have 


(Lo + 5lo)(yo + By)? = (11 + 8h) x? + 2(12 + l2)« + (1g + 5/8) (7) 
to the first order, and hence 


(Io +8Io)(co+ 8c)? = Ty +8; + 2(I2+8l2)/« + (1g + 8lg)/x2. (8) 


4- Conclusions 


The use of these results would be that it would be possible to examine the effects 
of small changes of p, A, ~ over non-overlapping intervals of z, and hence to use 
observed group-velocities to identify more closely the intervals where such changes 


are needed. They do for surface waves what Bullen’s transformation (1960) does 
for bodily waves. 
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The Nature of the Mohorovici¢ Discontinuity 


E. C. Bullard and D. T. Griggs 


(Received 1961 August 4) 


Summary 


It has been suggested that the Mohorovici¢ discontinuity represents 
a change in crystal structure and not a change in chemical composition. 
Three objections to this view are discussed and it is concluded that it is 
unlikely to be correct. ‘The objections are (1) that the stability relations 
are such that the continental and oceanic discontinuities cannot both 
have the heavier material below them; (2) that the small variation of 
the depth of the discontinuity with heat flow cannot be accounted for 
and (3) that the transition should not be sharp in a rock which is not 
composed of a single well defined chemical substance. 


1. Introduction 


It has often been supposed that the Mohorovitié discontinuity, or Moho, 
represents not the boundary between two different materials, but a change in 
crystal structure induced by increase of pressure with depth (Fermor 1914, 1938, 
Holmes 1927, Lovering 1958, Kennnedy 1959). This view is attractive since it 
suggests a mechanism that can cause widespread vertical movements. It is also 
plausible in that changes in crystal structure are observed in many silicates at 
about the pressures and temperatures occurring at the Moho. Unfortunately 
there are three serious objections, which it is the purpose of this paper to discuss. 

The first objection concerns the comparison between the variation of tempera- 
ture with depth in the Earth and the variation of the transition temperature with 
pressure. It is shown in Section 2 that the oceanic and continental transitions can 
not both result in the high pressure form lying beneath the Moho. 

The second objection (Section 3) is that the variations of Moho depth from 
place to place are less than would be expected from the variations in heat flow. 

The third objection (Section 4) is that a Moho produced by a phase change 
would be spread over a range of depth. 


2. The pressure-temperature relations 


‘The mean heat flow in the Atlantic and Pacific oceans away from the mid-ocean 
ridges has been estimated at 1-08 + 0°05 ycal/cm? s (Bullard 1961). This may be a 
little low since unpublished measurements by Von Herzen show higher values 
over a considerable area off the coast of California. On land the mean usually 
quoted is 1-2 ycal/cm? s, this may also be a little low since several recent measure- 
ments have been above it. We will assume 1-2 .cal/cm? s as a mean for both oceans 
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and continents. Beneath a continent the pressure and temperature have been 
computed on the assumption that the depth of the Moho is 35 km and that the 
material above it has a thermal conductivity of 0-006 cal/ems °C. Two extreme 
assumptions are considered, (1) that all the heat comes from below the Moho and 
(2) that it is all produced by heat sources uniformly distributed through the crust. 
The results are given in Table 1 which shows that the temperature at the Moho is 
likely to be between 350 and 7oo°C. As a rough approximation we take 
500 + 150 °C, this is equivalent to assuming that 57 per cent of the heat is produced 
in the crust. The temperature, 7, at any depth, D, in the crust is then given by 


T = 20D—0:'163D*". (1) 


Table 1 


Thick- Thermal Temperature 

Depth ness Density condy. Pressure Gradient (1) (2) 

km km g/cm® cal/ems°C kg/cm® °C/km °C “Cc 

Continent ° ° 
35 2°8 0006 


20 ° ° 
9800 86 700 350 
° == 


500 ° 
Sed. 


600 26 
Hard Rock 


1860 


The relation between temperature and pressure, p (in kg/cm?), is 


T = 0°:0714p— 2°08 x 10~6p?; 


this gives the line Continent in Figure 1, C representing the Moho. Under the 
ocean, above the Moho and below the sediments, 


T = 16°8+0°0714); 
this gives the line Ocean in Figure 1 where O represents the Moho. 

If the continental and oceanic Mohos are caused by the same transition, the 
temperature—pressure relation for the transition must be represented by some line 
through C and O. A straight line is shown in Figure 1 ; this line cuts the continental 
pressure temperature relation so that the denser material is stable below the Moho 
as it should be; under the ocean, however, the denser material is the stable phase 
at the ocean floor and the transition to the less dense form occurs at O. Provided 
that the continental and oceanic heat flows are approximately equal, this result 
does not depend critically on the numerical constants used nor on the transi- 
tion curve being a straight line. It is a consequence of the topology of Figure 1 
and requires only that the continental and oceanic temperature curves intersect 
above the continental Moho. It can be shown that this will happen if the con- 
tinental and oceanic heat flows are approximately equal and if more than 5 per 
cent of the continental heat flow is generated in the crust. 

If the transition line were strongly curved, as in Figure 2, it would intersect the 
oceanic pressure-temperature line in two points. The denser form of the material 
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would then be stable below the Moho, but would change again to the lighter form 
at a greater depth, O’. There is no evidence for such an ‘‘anti-Moho” and it does 
not seem likely that it exists. ‘The only escape from the conclusion that the contin- 
ental and oceanic Mohos cannot represent the same change in crystal structure 
appears to be in the not very plausible assumption that the measured oceanic 
heat flows are greatly disturbed by the intrusion of lava at shallow depths. 


Temperature (°C) 
200 400 
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Fic. 1.—The pressure-temperature relations in the crust and the 
transition line. 


3. The variation of Moho depth with heat flow 


The observed heat flows in non-volcanic districts on the continents vary from 
about 0°6 to 3-0 yxcal/cm? s; there are considerable areas where the mean heat flows 
are as low as 0°8 and as high as 2-4 ycal/cm? s. There appear to be no measurements 
of Moho depth in areas where the heat flow is known to be exceptionally high or low. 
It is very desirable that such measurements be made, but in their absence we assume 
that the usual spread of the observed values from 30 to 45 km has in fact covered 
some places of anomalous heat flow. 

If it is assumed that the temperatures in a continental region with heat flow H 
can be obtained from those given in Figure 1 by multiplying by H/1-2, the depth 
of the Moho can be found for any heat flow. For heat flows less than 1-05 cal/cm? s 
the temperature—pressure curve does not intersect the transition line of Figure 1 
and there is no solution. For larger heat flows the Moho depths are: 


Heat flow I'l 1°2 1°6 2°0 2°4 pcal/cm? s 
Moho depth 25 35 58 (71) (80) km 


For heat flows above 1-8 the solutions give temperatures decreasing with depth and 
are therefore presumably unacceptable. In an effort to widen the range of heat 
flows which can be accommodated, we drop the restriction that all temperature- 
depth curves in the crust must be similar and assume only that some part of the 
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heat is derived from radioactivity uniformly spread through the crust and that the 


rest comes from below the Moho. ‘The variation of temperature with pressure is 
then 


T = Bp—Cp. 


The coefficient B is fixed by the surface heat flow, Ho, and the density, p, and 
thermal conductivity, k, of crustal rocks, 


B = Ho/pk. 
The coefficient C then determines how great a proportion of the heat comes from 


below the Moho, it must be positive and not so large that the temperature at the 
Moho decreases downwards. 
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Fic. 2.—Continental and oceanic Mohos with strongly curved transition 
line. 


Families of temperature pressure curves for Hg = 0-8 and 2-4 are given in 
Figure 3 assuming h to be 0-006 ycal/cm?s° C and p to be 2-8 g/cm. The transi- 
tion curve must go near the point O, whose position is not seriously in doubt and 
must be cut by one curve from each family in the pressure range appropriate to the 
continental Moho (8 400 to 12 600 kg/cm? for depths between 30 and 45 km). The 
curves must have a positive slope at the points of intersection. From Figure 3 
it is seen that no such line exists even if moderate departures from straightness are 
permitted. 

The only way to avoid this difficulty appears to be to assume that in areas of 
high heat flow the sources of heat are at shallow depths in the crust and that in 
areas of low heat flow most of the heat comes from below the crust. These assump- 
tions would be quite arbitrary and do not seem particularly plausible when applied 
to such areas as Tasmania or Hungary where there are high heat flows and no signs 


of recent vulcanism. Measurements of Moho depths in these places are particularly 
desirable. 
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4. The sharpness of the transition 


The phase rule shows that if a pure substance exists in two crystal forms, the 
two forms can be in equilibrium along a line in the pressure temperature diagram 
as in Figure 1. That is, at a fixed temperature the transition takes place at a definite 
pressure. If the substance is capable of a continuous range of composition, as, for 
example, is the isomorphous series forsterite-fayalite, the two forms can be in 
equilibrium over an area in the pressure-temperature diagram. There is then no 
definite transition pressure, even if the temperature is held constant. The situa- 
tion is analogous to the melting of water. Pure water at atmospheric pressure 
freezes sharply at o °C, but if a salt solution is cooled, ice forms and remains in 
equilibrium with the brine over a range of temperature. 
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Fic. 3.—Temperature—pressure curves for heat flows of 2-4 (full lines) and 
0°8 peal/cm?s (dashed lines). O is the point representing the oceanic 
Moho from Figure 1. 


Since silicate minerals are nearly all members of isomorphous series, and are 
not compounds with a fixed composition, their structural transitions as the pres- 
sure is raised should not occur sharply, but should be spread over a range of 
pressure. Thus if the Moho represents a change in crystal structure it would not 
occur at a definite depth, but would be spread over a range, probably of many 
kilometres. Several workers claim to have observed reflections of seismic waves 
from the Moho (Steinhort & Meyer 1961), if these are genuine the interface cannot 
be spread over more than 2 or 3 kilometres. A more critical examination of these 
reflections is very desirable both on land and at sea. If they are genuine it is diffi- 
cult to believe that the Moho represents a phase change. 


5- Conclusions 


The state of the material at inaccessible depths within the Earth is necessarily 
subject to great uncertainty, and a wide range of assumptions can be made without 
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contradicting any observation. It does seem, however, that the most natural inter- 
pretation of the Moho as a phase change is untenable and that it is more likely 
that it represents a change in chemical composition. It is possible that some 
unexpected complication may vitiate the arguments of this paper, for example 
it might be that the material beneath the continental Moho is different from that 
beneath the oceanic Moho. Such doubts cannot be profitably discussed till we have 
more information, which should be provided in the foreseeable future by the drilling 
of oceanic Moholes. 

Just as this paper was being sent to press we received an MS of a paper by Dr G. 


Wetherill which is to be published in the J. Geophys. Res. and which contains many 
of the same ideas. 


Dept. of Geodesy & Geophysics, Institute of Geophysics, 
Madingley Rise, Los Angeles Laborataries, 
Madingley Road, Los Angeles, 
Cambridge: Cal., U.S.A. 
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Deformation of the Japanese Islands as inferred 
from Rock Magnetism 


N. Kawai,* H. Itot and S. Kume*t 
(Received 1960 December 21) 


Summary 


Results are presented of the directions of magnetization for twenty- 
eight sets of rocks collected from different parts of Japan. The directions 
for rocks of Pre-Tertiary age from south-west Japan differ systematically 
from those from north-east Japan. It is suggested that the difference is 
due to the deformation of the Japanese Islands in the late Mesozoic or 
early ‘Tertiary. Because of the deformation, it is not possible to infer 


the accurate position of the geomagnetic dipole with respect to Japan 
for Pre-'Tertiary times. 


1. Introduction 


Oriented rock samples have been collected over a wide area of the Japanese 
Islands of Hokkaido, Honshu, Shikoku and Kyushu. The collecting sites are 





@ Pre-Tertiary rocks with stable magnetization 
O Rocks with unsteble magnetization or uncertian age 


@ Post-Tertiary rocks with stable magnetization 








135° 





Fic. 1.—The location of sites where rocks were collected. ‘The numbered 
sites are for rocks of Pre-Tertiary age and the numbers correspond to 
those given in Table 1 and Figures 2 and 3. 
* Geological and Mineralogical Institute, Kyoto University, Japan. 
+ Institute of Earth Science, Kyoto University, Japan. 
} Physics Dept., King’s College, University of Durham, England. 
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shown in Figure 1. The rocks range in age from late Palaeozoic to Quaternary and 
include a variety of rock types. 


2. Directions of natural remanent magnetization 


Laboratory measurements of natural remanent magnetization (n.r.m.) are 
summarized in Tables 1 and 2. ‘The method of measurement and accuracy to be 
expected for individual measurements have been given by Kawai (1951). For the 
interpretation given later, it is convenient to group the results from north-east and 
south-west Japan separately. The directions of n.r.m. for different sites are shown 
on Schmidt’s equal area projections in Figures 2 to 5. Figures 2 and 3 show the 
directions for Pre-Tertiary rocks from north-east and south-west Japan respectively 


and Figures 4 and 5 those for Tertiary and Quaternary rocks from north-east and 
south-west Japan. 


Table 2 
N.R.M. of Tertiary and Quaternary Rocks 


Locality Rock type D I S 


From north-east parts of Japan 
Asama Andesite N 5°E +37° 17 
(Nagano) 
. Asama Andesite N 10E +49 
(Nagano) 
Showashinzan Basalt N 10E +57 
(Hokkaido) 
Arafune Porphyrite N 24E +49 
(Nagano) 
Towada Andesite N o +56 
(Aomori) 


From south-west parts of Japan 
Tomuro Basalt N 1E +48 
(Ishikawa) 
. Tomoru Andesite N 7W +50 
(Ishikawa) 
Shigisan Andesite N 4W +59 
(Nara) 
Sakurajima Andesite N 7E +45 
(Kagoshima) 
. Tojimbo Andesite +54 
(Fukui) 
Sugiyama Andesite N 39E +49 
(Fukui) 
Ibaragi Tuff 3 6W —51 
(Osaka) 
Genbudo Basalt S 13W —35 
(Hyogo) 


Mean declination of n.r.m. 
Mean inclination of n.r.m. (deg). 
Number of rock samples. 


Radius of 95 per cent circle of confidence (deg) (Fisher 1953). 
Fisher’s precision factor. 


Consideration of all the data for rocks likely to possess a stable magnetization 
and whose ages are known reveals three interesting facts. First, all rocks of Tertiary 
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and Quaternary age have directions of magnetization nearly parallel or anti- 
parallel to the present dipole field. These are shown in Figures 6c and 6d where all 
the results are plotted as normal directions. Secondly, the mean directions of 
magnetization of Pre-Tertiary (Mesozoic and late Palaeozoic) rocks from the north- 
eastern part of Japan are found to be north-westerly (Figure 6a) and thirdly, the 
mean direction for rocks of late Palaeozoic and early Mesozoic age from the south- 
west part of Japan are found to be north-easterly (Figure 6b). ‘The mean directions 
of magnetization are summarized in Figure 6e. It is seen that whilst the mean 
directions of magnetization of rocks of ‘Tertiary and Quaternary ages from both 
north-east and south-west Japan are close to the axis of the present dipole field, 


those of pre-Tertiary age from north-east Japan differ in declination from those in 
south-west Japan by about 40°. 











eee 


Fic. 7.—Map of the Japanese Islands. AA’ and BB’ are the axes of the 

south-west and north-east parts of Japan and the arrows represent the 

mean directions of n.r.m. of Pre-Tertiary rocks sampled from these two 
regions. 





3- Interpretation of results 


It is suggested that the difference in declination of 40° for the n.r.m. of rocks 
collected from south-west and north-east Japan is due to the deformation of the 
Japanese Islands. Figure 7 shows the position of a zone of plastic deformation 
known in Japan as the “fossa magna”. There are striking differences in geological 


structure on either side of the “fossa magna” and earthquakes and volcanism are 
associated with this zone. 








1. Miyemori 2. Tono 3.Tono 
(serpentine) (gabbro) (serpentine) 





4. Hitokabe 5. Yamaya 6.Takata 
(granite) (granite) (grano-diorite) 

7. Kinkazan 8. Tono 9. Nakoso 
(granite) (gabbro) (gabbro) 





Fic. 2.—Directions of n.r.m. of Pre-Tertiary (Mesozoic) rocks collected 
from the north-eastern part of Japan. Schmidt equal-area projections 
are used throughout and points are plotted on the lower hemisphere. 
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tions of n.r.m. of Pre-Tertiary (late Palaeozoic and Fic. 4.—Directions of n.r.m. of Tertiary and Quaternary rocks colle 
cks collected from the south-western part of Japan. from the north-eastern part of Japan. 
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Fic. 5.—Directions of n.r.m. of Tertiary and Quaternary rocks collected Fic. 6.—Directions of 1 
from the south-western part of Japan. east and south-west pa 
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1e present dipole field for Japan. 
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Fic. 8.—Schematic diagram of the bending of the Japanese Islands 

in relation to the hypothesis of polar wandering. (A) Polar wander path 

as deduced for North America (B) Polar wander curve as deduced for 
Europe and (C) polar wander curve as deduced for Japan. 








Fic. 9.—A possible interpretation of the direction of magnetization in 
Pre-Tertiary times by correcting for the bending of the Japanese 
Islands. 
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Tamaki (1955) in an analysis of earthquake data has shown the existence of 
high density material near the Izu Peninsula, that is, near the “fossa magna’’. 
This high density sector may be considered to be a mass of basic rocks intruded 
into the fracture zones which resulted from the bending of the Japanese Islands 
along this line. Further, Ehara (1961) found gabbro and granodiorite associated 
with the andesitic magma extruded at Oshima Island confirming the existence of 
high density material. The region is also associated with positive gravity anomalies 
(Matsuyama 1936, and Kumagai 1940 and 1953). 

Two lines, AA’ and BB’ in Figure 7 represent the mean axes of south-west and 
north-east Japan respectively. The angle between these lines (@) is approximately 
equal to the difference in the mean directions of magnetization for rocks of pre- 
Tertiary age collected from these two parts of Japan. It is suggested that the 
difference in direction of magnetization for rocks from north-east and south-west 
of Japan is due to the bending of the Japanese Islands along the line of the “‘fossa 
magna” and in agreement with the geological evidence. Further, it is possible to 
infer from the rock magnetism data that the bending took place in early Tertiary 
times. 


4- Implications for studies of polar wandering 


In view of this rotational movement of the different parts of the Japanese 
Islands, it is very difficult to infer the position of the geomagnetic pole with respect 
to Japan for pre-Tertiary times. If the amount of polar wander between the 
Mesozoic and early ‘Tertiary is assumed to be small, it is possible to infer that the 
amount of relative rotation of north-east and south-west Japan is of the order of 


° 


40 


The polar wander curve for Japan as deduced by Nagata & others (1959) 
differs greatly from those deduced for Europe and North America by Creer & 
others (1954) and by Runcorn (1959) as shown in Figure 8. The discrepancy 
might be considerably reduced if the bending of the Japanese Islands is taken 
into consideration. Figure 9 shows one possible interpretation of removing the 
effect of the bending of the Japanese Islands and the difficulty of defining the 
polar wander curve for Japan in pre-Tertiary times. 


5- Acknowledgments 


The authors wish to thank Professor N. Kumagai of Kyoto University for 
his continuous encouragement throughout this study. Their thanks are also due 
to Dr R. W. Girdler of Durham University for his kind suggestions and 
correction of English. 


Physics Dept., King’s College, 
University of Durham, 
Newcastle-upon- Tyne: 
1960 December. 





N. Kawai, H. Ito and S. Kume 


References 


Creer, K. M., Irving E., & Runcorn, S. K., 1954. J. Geomag. & Geoelect., 6, 163. 

Ehara, S., 1961. Geogr. Rev., Jap., 34 1. 

Fisher, Sir Ronald, 1953. Proc. Roy. Soc. A, 217, 295. 

Kawai, N., 1951. J. Geophys. Res., 56, 73. 

Kumagai, N., 1940. Jap. J. Astron. Geophys. 17, 477. 

Kumagai, N., 1953. Rep. Jap. Geod. Commis., No. 4. 

Matsuyama, M., 1936. Proc. Imp. Acad. Jap., 12, 93. 

Nagata, T., Akimoto, S., Shimizu, Y., Kobayashi, K. and Kuno, H., 1959. Proc. 
Jap. Acad., 35, 378. 

Runcorn, S. K., 1959. Science, 129, 1002. 

Tamaki, I., 1955. Zisin, 7, 226. 





N. Kawai, H. Ito and S. Kume 


References 


Creer, K. M., Irving E., & Runcorn, S. K., 1954. J. Geomag. & Geoelect., 6, 163. 

Ehara, S., 1961. Geogr. Rev., Jap., 34, 1. 

Fisher, Sir Ronald, 1953. Proc. Roy. Soc. A, 217, 295. 

Kawai, N., 1951. J. Geophys. Res., 56, 73. 

Kumagai, N., 1940. Jap. J. Astron. Geophys. 17, 477. 

Kumagai, N., 1953. Rep. Jap. Geod. Commis., No. 4. 

Matsuyama, M., 1936. Proc. Imp. Acad. Jap., 12, 93. 

Nagata, T., Akimoto, S., Shimizu, Y., Kobayashi, K. and Kuno, H., 1959. Proc. 
Jap. Acad., 35, 378. 

Runcorn, S. K., 1959. Science, 129, 1002. 

Tamaki, I., 1955. Zisin, 7, 226. 





